06_Analyses_on_Hospital_Discharge

# Load packages

# Ensures the package "pacman" is installed
if (!require("pacman")) install.packages("pacman")

# p_load() installs (if necessary) and loads the following packages

pacman::p_load(plyr, # for the Varying baseline risk calculations
               magrittr, # same as above
               rio, # to import/export files
               tidyverse, # Data wrangling and plotting
               kableExtra, # To create tables
               gt, # To create tables
               flextable, # To create tables
               here, # To find file paths within R project
               metafor, # To create priors from meta-analysis
               janitor, # To change column names
               ggdist, # To plot distributions
               patchwork, # To arrange plots
               ggsci, # Color palettes
               Cairo # To export PDF files
               )

# Load data on respiratory support subgroups

d1 = import((here("final_analyses", "output", "data", # file path
                   "respiratory-data.csv")),           # file name
             header = T)

# Load data on corticosteroids subgroups

d2 = import(here("final_analyses", "data", # file path
                            "corticosteroids_subgroup_extracted-data.xlsx")) 

d2 = clean_names(d2)

In this document, we will compile most of the analyses regarding hospital discharge. First, we will derive the risk difference from odds ratio. Next, we will perform the main analyses and sensitivity analyses. We will be brief, so for complete explanation on the rationale of the analyses, please check the files on the mortality outcome.

In summary, this document shows what all other documents has shown, but now the spotlight is on the hospital discharge outcome.

Of note, for this outcome, positive log-odds ratio, odds ratio greater than 1, and negative risk differences mean benefit.

Calculating posterior distributions

Priors and RECOVERY

Use of corticosteroids

Not using

First, let’s calculate and visualize the log-odds ratio from the RECOVERY trial.

logOR_recovery_not_using =
  rma(
    # Outcome
    measure = "OR", # Log-odds ratio
    
    # Tocilizumab group
    ai = events_toci,
    n1i = total_toci,

    # Control group
    ci = events_control,
    n2i = total_control,
    
    data = d2 %>% filter(outcome == "discharge",
                         corticosteroids == 0,
                         trial == "RECOVERY"),
    
     method = "REML",

    slab = paste(trial)
  )

forest(logOR_recovery_not_using)

Now, let’s calculate and visualize the log-odds ratio from the other trials. The pooled result from this random-effect meta-analysis will be used as the estimate for the distribution of prior regarding this subgroup.

# Simple oxygen

logOR_prior_not_using =
  rma(
    # Outcome
    measure = "OR", # Log-odds ratio
    
    # Tocilizumab group
    ai = events_toci,
    n1i = total_toci,

    # Control group
    ci = events_control,
    n2i = total_control,
    
    data = d2 %>% filter(outcome == "discharge",
                         corticosteroids == 0,
                         trial != "RECOVERY"),
    
     method = "REML",

    slab = paste(trial)
  )

forest(logOR_prior_not_using)

Using

Let’s calculate and visualize the log-odds ratio from the RECOVERY trial.

logOR_recovery_using =
  rma(
    # Outcome
    measure = "OR", # Log-odds ratio
    
    # Tocilizumab group
    ai = events_toci,
    n1i = total_toci,

    # Control group
    ci = events_control,
    n2i = total_control,
    
    data = d2 %>% filter(outcome == "discharge",
                         corticosteroids == 1,
                         trial == "RECOVERY"),
    
     method = "REML",

    slab = paste(trial)
  )

forest(logOR_recovery_using)

Now, let’s calculate and visualize the log-odds ratio from the other trials. The pooled result from this random-effect meta-analysis will be used as the estimate for the distribution of prior regarding this subgroup.

# Simple oxygen

logOR_prior_using =
  rma(
    # Outcome
    measure = "OR", # Log-odds ratio
    
    # Tocilizumab group
    ai = events_toci,
    n1i = total_toci,

    # Control group
    ci = events_control,
    n2i = total_control,
    
    data = d2 %>% filter(outcome == "discharge",
                         corticosteroids == 1,
                         trial != "RECOVERY"),
    
     method = "REML",

    slab = paste(trial)
  )

forest(logOR_prior_using)

Respiratory support subgroups

Simple oxygen only

Let’s calculate and visualize the log-odds ratio from the RECOVERY trial.

logOR_recovery_simple =
  rma(
    # Outcome
    measure = "OR", # Log-odds ratio
    
    # Tocilizumab group
    ai = events_toci,
    n1i = total_toci,

    # Control group
    ci = events_control,
    n2i = total_control,
    
    data = d1 %>% filter(outcome == "discharge",
                         oxygen == "simple oxygen",
                         trial == "RECOVERY"),
    
     method = "REML",

    slab = paste(trial)
  )

forest(logOR_recovery_simple)

Now, let’s calculate and visualize the log-odds ratio from the other trials. The pooled result from this random-effect meta-analysis will be used as the estimate for the distribution of prior regarding this subgroup.

logOR_prior_simple =
  rma(
    # Outcome
    measure = "OR", # Log-odds ratio
    
    # Tocilizumab group
    ai = events_toci,
    n1i = total_toci,

    # Control group
    ci = events_control,
    n2i = total_control,
    
    data = d1 %>% filter(outcome == "discharge",
                         oxygen == "simple oxygen",
                         trial != "RECOVERY"),
    
     method = "REML",

    slab = paste(trial)
  )

forest(logOR_prior_simple)

Non-invasive ventilation

Let’s calculate and visualize the log-odds ratio from the RECOVERY trial.

logOR_recovery_noninvasive =
  rma(
    # Outcome
    measure = "OR", # Log-odds ratio
    
    # Tocilizumab group
    ai = events_toci,
    n1i = total_toci,

    # Control group
    ci = events_control,
    n2i = total_control,
    
    data = d1 %>% filter(outcome == "discharge",
                         oxygen == "non-invasive ventilation",
                         trial == "RECOVERY"),
    
     method = "REML",

    slab = paste(trial)
  )

forest(logOR_recovery_noninvasive)

Now, let’s calculate and visualize the log-odds ratio from the other trials. The pooled result from this random-effect meta-analysis will be used as the estimate for the distribution of prior regarding this subgroup.

logOR_prior_noninvasive =
  rma(
    # Outcome
    measure = "OR", # Log-odds ratio
    
    # Tocilizumab group
    ai = events_toci,
    n1i = total_toci,

    # Control group
    ci = events_control,
    n2i = total_control,
    
    data = d1 %>% filter(outcome == "discharge",
                         oxygen == "non-invasive ventilation",
                         trial != "RECOVERY"),
    
     method = "REML",

    slab = paste(trial)
  )

forest(logOR_prior_noninvasive)

Invasive mechanical ventilation

Let’s calculate and visualize the log-odds ratio from the RECOVERY trial.

logOR_recovery_invasive =
  rma(
    # Outcome
    measure = "OR", # Log-odds ratio
    
    # Tocilizumab group
    ai = events_toci,
    n1i = total_toci,

    # Control group
    ci = events_control,
    n2i = total_control,
    
    data = d1 %>% filter(outcome == "discharge",
                         oxygen == "invasive ventilation",
                         trial == "RECOVERY"),
    
     method = "REML",

    slab = paste(trial)
  )

forest(logOR_recovery_invasive)

Now, let’s calculate and visualize the log-odds ratio from the other trials. The pooled result from this random-effect meta-analysis will be used as the estimate for the distribution of prior regarding this subgroup.

logOR_prior_invasive =
  rma(
    # Outcome
    measure = "OR", # Log-odds ratio
    
    # Tocilizumab group
    ai = events_toci,
    n1i = total_toci,

    # Control group
    ci = events_control,
    n2i = total_control,
    
    data = d1 %>% filter(outcome == "discharge",
                         oxygen == "invasive ventilation",
                         trial != "RECOVERY"),
    
     method = "REML",

    slab = paste(trial)
  )

forest(logOR_prior_invasive)

Posterior Distribution

## Load custom functions

# To generate normal conjugate posterior probabilities

source(here("final_analyses", "script", "functions", "analyses", # file path
            "conjugate_normal_posterior.R"))                     # file name

# To generate posterior probabilities with an evidence-based prior
# (it requires "conjugate_normal_posterior.R" )

source(here("final_analyses", "script", "functions", "analyses", # file path
            "normal_approximation_logOR.R"))                     # file name

# To generate posterior probabilities with multiple priors
# (it requires "conjugate_normal_posterior.R" )

source(here("final_analyses", "script", "functions", "analyses", # file path
            "logOR_normal_approximation_multiple_priors.R"))     # file name

# To compile all posteriors into a tibble

source(here("final_analyses", "script", "functions", "analyses", # file path
            "tibble_all_posteriors_logOR.R"))                      # file name

Note that positive log(OR) means benefit for the hospital discharge outcome. Thus, in this document, the optimistic prior has a positive mean and the pessimistic the opposite.

## Use of corticosteroids

# Not using
posteriors_not_using_discharge =
  logOR_normal_approximation_multiple_priors(
  # Data object with evidence-based posterior from normal_approximation()
  posterior_not_using_discharge, 
  "discharge", # Outcome (within quotes)
  posteriors_not_using_discharge # Name for the final output with data+prior+posterior
  )

## Post mean =  -0.01606234 , Post Var =  0.0222785 , Post SD =  0.1492598 
## Post mean =  0.09201482 , Post Var =  0.01767745 , Post SD =  0.1329566 
## Post mean =  -0.01274508 , Post Var =  0.01767745 , Post SD =  0.1329566 
## Post mean =  0.03333988 , Post Var =  0.01767745 , Post SD =  0.1329566 
## Post mean =  -0.05883004 , Post Var =  0.01767745 , Post SD =  0.1329566

# Using
posteriors_using_discharge =
  logOR_normal_approximation_multiple_priors(
  # Data object with evidence-based posterior from normal_approximation()
  posterior_using_discharge, 
  "discharge", # Outcome (within quotes)
  posteriors_using_discharge # Name for the final output with data+prior+posterior
  )

## Post mean =  0.3479437 , Post Var =  0.004814994 , Post SD =  0.06939016 
## Post mean =  0.3180544 , Post Var =  0.004571678 , Post SD =  0.06761418 
## Post mean =  0.330361 , Post Var =  0.004571678 , Post SD =  0.06761418 
## Post mean =  0.3416372 , Post Var =  0.004571678 , Post SD =  0.06761418 
## Post mean =  0.3190848 , Post Var =  0.004571678 , Post SD =  0.06761418

## Respiratory support

# Simple oxygen only
posteriors_simple_discharge =
  logOR_normal_approximation_multiple_priors(
  # Data object with evidence-based posterior from normal_approximation()
  posterior_simple_discharge, 
  "discharge", # Outcome (within quotes)
  posteriors_simple_discharge # Name for the final output with data+prior+posterior
  )

## Post mean =  0.3179829 , Post Var =  0.0105669 , Post SD =  0.1027954 
## Post mean =  0.331199 , Post Var =  0.009596272 , Post SD =  0.09796056 
## Post mean =  0.2887745 , Post Var =  0.009596272 , Post SD =  0.09796056 
## Post mean =  0.3092716 , Post Var =  0.009596272 , Post SD =  0.09796056 
## Post mean =  0.2682773 , Post Var =  0.009596272 , Post SD =  0.09796056

# Non-invasive ventilation
posteriors_noninvasive_discharge =
  logOR_normal_approximation_multiple_priors(
  # Data object with evidence-based posterior from normal_approximation()
  posterior_noninvasive_discharge, 
  "discharge", # Outcome (within quotes)
  posteriors_noninvasive_discharge # Name for the final output with data+prior+posterior
  )

## Post mean =  0.2913939 , Post Var =  0.00962873 , Post SD =  0.09812609 
## Post mean =  0.2883111 , Post Var =  0.008837633 , Post SD =  0.09400869 
## Post mean =  0.267453 , Post Var =  0.008837633 , Post SD =  0.09400869 
## Post mean =  0.2857867 , Post Var =  0.008837633 , Post SD =  0.09400869 
## Post mean =  0.2491193 , Post Var =  0.008837633 , Post SD =  0.09400869

# Invasive mechanical ventilation
posteriors_invasive_discharge =
  logOR_normal_approximation_multiple_priors(
  # Data object with evidence-based posterior from normal_approximation()
  posterior_invasive_discharge, 
  "discharge", # Outcome (within quotes)
  posteriors_invasive_discharge # Name for the final output with data+prior+posterior
  )

## Post mean =  0.2352049 , Post Var =  0.04918534 , Post SD =  0.2217777 
## Post mean =  0.2702109 , Post Var =  0.03654003 , Post SD =  0.1911545 
## Post mean =  0.1747349 , Post Var =  0.03654003 , Post SD =  0.1911545 
## Post mean =  0.2321048 , Post Var =  0.03654003 , Post SD =  0.1911545 
## Post mean =  0.117365 , Post Var =  0.03654003 , Post SD =  0.1911545

## Use of corticosteroids

# Not using
samples_logOR_posteriors_not_using_discharge = 
  tibble_all_posteriors_logOR(posteriors_not_using_discharge)

# Using
samples_logOR_posteriors_using_discharge = 
  tibble_all_posteriors_logOR(posteriors_using_discharge)

## Respiratory support

# Simple oxygen only
samples_logOR_posteriors_simple_discharge = 
  tibble_all_posteriors_logOR(posteriors_simple_discharge)

# Non-invasive ventilation
samples_logOR_posteriors_noninvasive_discharge = 
  tibble_all_posteriors_logOR(posteriors_noninvasive_discharge)

# Invasive mechanical ventilation
samples_logOR_posteriors_invasive_discharge = 
  tibble_all_posteriors_logOR(posteriors_invasive_discharge)

## Use of corticosteroids

# Not using
samples_posteriors_not_using_discharge = 
  samples_logOR_posteriors_not_using_discharge %>% 
  pivot_longer(
    cols = c('Non-informative':'Pessimistic'),
    names_to = "Underlying_Prior", values_to = "log(OR)"
  )

# Using
samples_posteriors_using_discharge = 
  samples_logOR_posteriors_using_discharge %>% 
  pivot_longer(
    cols = c('Non-informative':'Pessimistic'),
    names_to = "Underlying_Prior", values_to = "log(OR)"
  )

## Respiratory support

# Simple oxygen only
samples_posteriors_simple_discharge = 
  samples_logOR_posteriors_simple_discharge %>% 
  pivot_longer(
    cols = c('Non-informative':'Pessimistic'),
    names_to = "Underlying_Prior", values_to = "log(OR)"
  )

# Non-invasive ventilation
samples_posteriors_noninvasive_discharge = 
  samples_logOR_posteriors_noninvasive_discharge %>% 
  pivot_longer(
    cols = c('Non-informative':'Pessimistic'),
    names_to = "Underlying_Prior", values_to = "log(OR)"
  )

# Invasive mechanical ventilation
samples_posteriors_invasive_discharge = 
  samples_logOR_posteriors_invasive_discharge %>% 
  pivot_longer(
    cols = c('Non-informative':'Pessimistic'),
    names_to = "Underlying_Prior", values_to = "log(OR)"
  )

## Use of corticosteroids

# Not using
samples_posteriors_not_using_discharge = 
  samples_posteriors_not_using_discharge %>% 
  mutate(OR = exp(`log(OR)`))

# Using
samples_posteriors_using_discharge = 
  samples_posteriors_using_discharge %>% 
  mutate(OR = exp(`log(OR)`))

## Respiratory support

# Simple oxygen only
samples_posteriors_simple_discharge = 
  samples_posteriors_simple_discharge %>% 
  mutate(OR = exp(`log(OR)`))

# Non-invasive ventilation
samples_posteriors_noninvasive_discharge = 
  samples_posteriors_noninvasive_discharge %>% 
  mutate(OR = exp(`log(OR)`))

# Invasive mechanical ventilation
samples_posteriors_invasive_discharge = 
  samples_posteriors_invasive_discharge %>% 
  mutate(OR = exp(`log(OR)`))

Baseline risks in RECOVERY

Remember that higher risks of hospital discharge is a benefitial outcome.

Subgroup	Baseline Risk in Control Group (%)
Use of corticosteroids
Not using	45.8
Using	50.7
Respiratory support
Simple oxygen only	68.1
Non-invasive ventilation	41.8
Invasive mechanical	16.0

# BR = Baseline Risk
# OR = Odds Ratio

RD_fun = function(BR,OR)
  {
( BR*(BR-1)*(OR-1) ) / ( BR*OR - BR+1)
}

# Calculating the risk differences

### Use of corticosteroids

## Not using

# Define baseline risk
BR_not_using = d2 %>% 
    filter(trial == "RECOVERY",
           corticosteroids == 0, # Not using
           outcome == "discharge") %>% 
      summarise(n = events_control/total_control) %>%
      pull()

samples_posteriors_not_using_discharge = 
  samples_posteriors_not_using_discharge %>% 
  # Apply formula from Doi et al., 2020
  mutate(RD = RD_fun(BR_not_using,
                     OR))

## Using

# Define baseline risk
BR_using = d2 %>% 
    filter(trial == "RECOVERY",
           corticosteroids == 1, # Using
           outcome == "discharge") %>% 
      summarise(n = events_control/total_control) %>%
      pull()

samples_posteriors_using_discharge = 
  samples_posteriors_using_discharge %>% 
  # Apply formula from Doi et al., 2020
  mutate(RD = RD_fun(BR_using,
                     OR))

### Respiratory support

## Simple oxygen

# Define baseline risk
BR_simple = d1 %>% 
    filter(trial == "RECOVERY",
           oxygen == "simple oxygen",
           outcome == "discharge") %>% 
      summarise(n = events_control/total_control) %>% 
      pull()

samples_posteriors_simple_discharge = 
  samples_posteriors_simple_discharge %>% 
  # Apply formula from Doi et al., 2020
  mutate(RD = RD_fun(BR_simple,
                     OR))

## Non-invasive

# Define baseline risk
BR_noninvasive = d1 %>% 
    filter(trial == "RECOVERY",
           oxygen == "non-invasive ventilation",
           outcome == "discharge") %>% 
      summarise(n = events_control/total_control) %>% 
      pull()

samples_posteriors_noninvasive_discharge = 
  samples_posteriors_noninvasive_discharge %>% 
  # Apply formula from Doi et al., 2020
  mutate(RD = RD_fun(BR_noninvasive,
                     OR))

## Invasive mechanical ventilation

# Define baseline risk
BR_invasive = d1 %>% 
    filter(trial == "RECOVERY",
           oxygen == "invasive ventilation",
           outcome == "discharge") %>% 
      summarise(n = events_control/total_control) %>% 
      pull()

samples_posteriors_invasive_discharge = 
  samples_posteriors_invasive_discharge %>% 
  # Apply formula from Doi et al., 2020
  mutate(RD = RD_fun(BR_invasive,
                     OR))

Summary

To display some of our calculations, we will now plot the posterior distributions for log-odds ratio, odds ratio and risk difference for a single subgroup (not using corticosteroids).

We note that log-odds ratios greater than \(0\), odds-ratios greater than \(1\) and risk differences lower than \(0\) mean benefit for the hospital discharge outcome.

Analyses

We are going to replicate the same analyses done regarding the mortality outcome. Of course, now with the hospital discharge data.

One could load all data objects generated above without the need of running all code chunks again. THe purpose of the following code chunk is to facilitate to this end, by loading the data objects stored elsewhere.

## Load objects with mean + SD of multiple priors along with
## RECOVERY data and posterior distributions

#load(file = here("final_analyses", "output", "data", "analyses", # file path
#                 "data_prior_posteriors_mean_sd_discharge.RData"))   # file name

## Load posterior distributions' samples with logOR, OR and RD

# I uploaded these data objects to OSF since the .RData file is 56MB.
# Here I load all data files directly from OSF without having to download it
# to my/your PC

#samples = url("https://osf.io/7t32f/download?version=1")
#load(gzcon(samples))

Main figures and tables

# Load original data extraction table

d = import((here("final_analyses", "output", "data", # file path
                   "respiratory-data.csv")),           # file name
             header = T)
d = d %>% 
  as_tibble() %>% 
  select(!starts_with("X")) 

# Create different variables per subgroup

priors_simple_oxygen_raw =
  d %>%
  filter(
    trial != "RECOVERY",
    oxygen == "simple oxygen",
    outcome == "discharge"
  )

priors_noninvasive_raw =
  d %>%
  filter(
    trial != "RECOVERY",
    oxygen == "non-invasive ventilation",
    outcome == "discharge"
  )

priors_invasive_raw =
  d %>%
  filter(
    trial != "RECOVERY",
    oxygen == "invasive ventilation",
    outcome == "discharge"
  )

recovery_simple_oxygen = d %>%
  filter(trial == "RECOVERY",
         outcome == "discharge",
         oxygen == "simple oxygen")

recovery_noninvasive = d %>%
  filter(trial == "RECOVERY",
         outcome == "discharge",
         oxygen == "non-invasive ventilation")

recovery_invasive = d %>%
  filter(trial == "RECOVERY",
         outcome == "discharge",
         oxygen == "invasive ventilation")
  

# Load corticosteroids raw data

d_cortico = import(here("final_analyses", "data", # file path
                            "corticosteroids_subgroup_extracted-data.xlsx"))

d_cortico = clean_names(d_cortico)

priors_not_using =
  d_cortico %>%
  filter(
    trial != "RECOVERY",
    corticosteroids == 0,
    outcome == "discharge"
  )

priors_using =
  d_cortico %>%
  filter(
    trial != "RECOVERY",
    corticosteroids == 1,
    outcome == "discharge"
  )

recovery_not_using = d_cortico %>%
  filter(trial == "RECOVERY",
         outcome == "discharge",
         corticosteroids == 0)

recovery_using = d_cortico %>%
  filter(trial == "RECOVERY",
         outcome == "discharge",
         corticosteroids == 1)

# Set seed for reproducibility

set.seed = 123 

Overview in priors and RECOVERY

respiratory_table = 
  d %>% 
  # only mortality
  filter(outcome == "discharge") %>% 
  # change names to display them in the table
  mutate(oxygen = case_when(
    oxygen == "simple oxygen" ~ "Simple oxygen only",
    oxygen == "non-invasive ventilation" ~ "Non-invasive ventilation",
    oxygen == "invasive ventilation" ~ "Invasive mechanical ventilation",
  )) %>% 
  # Remove any row with NA, eg, the ones with "Pooled" in oxygen
  drop_na() %>%
  # Create columns with Risk (%) for both control and tocilizumab
  mutate(risk_control = 100*round(events_control/total_control, 2),
         risk_toci = 100*round(events_toci/total_toci,2)) %>% 
  # Order columns
  select(oxygen,
         trial,
         events_control, total_control, risk_control,
         events_toci, total_toci, risk_toci) %>% 
  # Rename to match the next table
  rename("Subgroup" = oxygen)

cortico_table = 
  d_cortico %>% 
  filter(outcome == "discharge") %>% 
  mutate(corticosteroids = case_when(
    corticosteroids == 0 ~ "Not using corticosteroids",
    corticosteroids == 1 ~ "Using corticosteroids"
  )) %>%
  mutate(risk_control = 100*round(events_control/total_control, 2),
         risk_toci = 100*round(events_toci/total_toci,2)) %>% 
  select(corticosteroids,
         trial,
         events_control, total_control, risk_control,
         events_toci, total_toci, risk_toci) %>% 
  rename("Subgroup" = corticosteroids)

# Put both tables togehter
both_tables = bind_rows(respiratory_table, cortico_table)

table1 = both_tables %>% 
  # Rename for final table
  rename("Study" = trial,
         "Events" = events_control,
         "Total" = total_control,
         "Risk (%)" = risk_control,
         "Events " = events_toci,
         "Total " = total_toci,
         "Risk (%) " = risk_toci) %>% 
  # Define order for final table
  mutate(Subgroup = factor(Subgroup,
                        levels = c("Not using corticosteroids",
                                   "Using corticosteroids",
                                   "Simple oxygen only",
                                   "Non-invasive ventilation",
                                   "Invasive mechanical ventilation")),
         Study = factor(Study,
                        levels = c("RECOVERY",
                                   "COVACTA",
                                   "REMAP-CAP",
                                   "CORIMUNO-19",
                                   "Salvarini"))) %>%
  # Re-arrange based on Study name
  arrange(Study) %>% 
  # Gropup to create row tabs with gt()
  group_by(Subgroup) %>% 
  # https://stackoverflow.com/questions/43832434/arrange-within-a-group-with-dplyr
  arrange(Subgroup, .by_group = TRUE) %>% 
  
  ### Transform into HTML table with gt()
  gt() %>% 
  # Centralize columns
  cols_align(
    align = "center",
    columns = everything()
  ) %>% 
  # Expand the size of first column
  cols_width(
    Study ~ px(250)
  ) %>% 
  # Add column tabs based on the column names defined above
  tab_spanner(label = "Control",
              columns = c("Events", "Total", "Risk (%)")) %>% 
  tab_spanner(label = "Tocilizumab",
              columns = c("Events ", "Total ", "Risk (%) "))

table1 

Study	Control			Tocilizumab
	Events	Total	Risk (%)	Events	Total	Risk (%)
Not using corticosteroids
RECOVERY	168	367	46	162	357	45
COVACTA	35	65	54	124	188	66
Using corticosteroids
RECOVERY	873	1721	51	987	1664	59
COVACTA	36	79	46	42	106	40
Simple oxygen only
RECOVERY	635	933	68	697	935	75
COVACTA	35	44	80	66	78	85
CORIMUNO-19	49	67	73	52	63	83
Non-invasive ventilation
RECOVERY	362	867	42	401	819	49
COVACTA	19	39	49	56	94	60
Salvarini	58	63	92	54	60	90
Invasive mechanical ventilation
RECOVERY	47	294	16	52	268	19
COVACTA	13	55	24	35	113	31

# table1 %>% 
#   gtsave(here("final_analyses", "output", "plots", "appendix", # File path
#                    "STable8_discharge_table1.png"))

## DEPRECATED

t = tibble(
  
  "Subgroup" = c(
    "Not using",
    "Using",
    
    "Simple oxygen only",
    "Non-invasive ventilation",
    "Invasive mechanical ventilation"
    
  ),
  
  ## Priors
  
  "Number of studies" = c(
    
    priors_not_using %>%
      count() %>% pull(),
    priors_using %>%
      count() %>% pull(),
    
    priors_simple_oxygen_raw %>%
      count() %>% pull(),
    priors_noninvasive_raw %>%
      count() %>% pull(),
    priors_invasive_raw %>%
      count() %>% pull()
    
  ),
  
  # Control arm
  
  "Events  " = c(
    
    priors_not_using %>%
    summarise(n = sum(events_control)) %>% pull(),
    priors_using %>%
      summarise(n = sum(events_control)) %>% pull(),
    
    priors_simple_oxygen_raw %>%
      summarise(n = sum(events_control)) %>% pull(),
    priors_noninvasive_raw %>%
      summarise(n = sum(events_control)) %>% pull(),
    priors_invasive_raw %>%
      summarise(n = sum(events_control)) %>% pull()

  ),
  
  "Total  " = c(
    
    priors_not_using %>%
      summarise(n = sum(total_control)) %>% pull(),
    priors_using %>%
      summarise(n = sum(total_control)) %>% pull(),
    
    priors_simple_oxygen_raw %>%
      summarise(n = sum(total_control)) %>% pull(),
    priors_noninvasive_raw %>%
      summarise(n = sum(total_control)) %>% pull(),
    priors_invasive_raw %>%
      summarise(n = sum(total_control)) %>% pull()

  ),
  
   "Risk (%)" = c(
    
    priors_not_using %>%
      summarise(n = round(100*sum(events_control)/sum(total_control),1)) %>%
      pull(),
    priors_using %>%
      summarise(n = round(100*sum(events_control)/sum(total_control),1)) %>%
      pull(),
    
    priors_simple_oxygen_raw %>%
      summarise(n = round(100*sum(events_control)/sum(total_control),1)) %>%
      pull(),
    priors_noninvasive_raw %>%
      summarise(n = round(100*sum(events_control)/sum(total_control),1)) %>%
      pull(),
    priors_invasive_raw %>%
      summarise(n = round(100*sum(events_control)/sum(total_control),1)) %>%
      pull()

  ),
  
  # Tocilizumab arm
  
  "Events   " = c(
    priors_not_using %>%
      summarise(n = sum(events_toci)) %>% pull(),
    priors_using %>%
      summarise(n = sum(events_toci)) %>% pull(),
    
    priors_simple_oxygen_raw %>%
      summarise(n = sum(events_toci)) %>% pull(),
    priors_noninvasive_raw %>%
      summarise(n = sum(events_toci)) %>% pull(),
    priors_invasive_raw %>%
      summarise(n = sum(events_toci)) %>% pull()
    
  ),
  
  "Total   " = c(
    priors_not_using %>%
      summarise(n = sum(total_toci)) %>% pull(),
    priors_using %>%
      summarise(n = sum(total_toci)) %>% pull(),
    
    priors_simple_oxygen_raw %>%
      summarise(n = sum(total_toci)) %>% pull(),
    priors_noninvasive_raw %>%
      summarise(n = sum(total_toci)) %>% pull(),
    priors_invasive_raw %>%
      summarise(n = sum(total_toci)) %>% pull()

  ),
  
  "Risk (%) " = c(
  
  priors_not_using %>%
    summarise(n = round(100*sum(events_toci)/sum(total_toci),1)) %>%
    pull(),
  priors_using %>%
    summarise(n = round(100*sum(events_toci)/sum(total_toci),1)) %>%
    pull(),
  
  priors_simple_oxygen_raw %>%
    summarise(n = round(100*sum(events_toci)/sum(total_toci),1)) %>%
    pull(),
  priors_noninvasive_raw %>%
    summarise(n = round(100*sum(events_toci)/sum(total_toci),1)) %>%
    pull(),
  priors_invasive_raw %>%
    summarise(n = round(100*sum(events_toci)/sum(total_toci),1)) %>%
    pull()
  
),
  
  ## RECOVERY
  
  # Control arm
  
  "Events" = c(
    recovery_not_using %>% 
      select(events_control) %>% pull(),
    recovery_using %>% 
      select(events_control) %>% pull(),
    
    recovery_simple_oxygen %>% 
      select(events_control) %>% pull(),
    recovery_noninvasive %>% 
      select(events_control) %>% pull(),
    recovery_invasive %>% 
      select(events_control) %>% pull()

  ),
  
  "Total" = c(
    
    recovery_not_using %>% 
      select(total_control) %>% pull(),
    recovery_using %>% 
      select(total_control) %>% pull(),
    
    recovery_simple_oxygen %>% 
      select(total_control) %>% pull(),
    recovery_noninvasive %>% 
      select(total_control) %>% pull(),
    recovery_invasive %>% 
      select(total_control) %>% pull()

  ),

  "Risk (%)  " = c(
    
    recovery_not_using %>% 
      summarise(n = round(100*sum(events_control)/sum(total_control),1)) %>%
      pull(),
    recovery_using %>% 
      summarise(n = round(100*sum(events_control)/sum(total_control),1)) %>%
      pull(),
    
    recovery_simple_oxygen %>% 
      summarise(n = round(100*sum(events_control)/sum(total_control),1)) %>%
      pull(),
    recovery_noninvasive %>% 
      summarise(n = round(100*sum(events_control)/sum(total_control),1)) %>%
      pull(),
    recovery_invasive %>% 
      summarise(n = round(100*sum(events_control)/sum(total_control),1)) %>%
      pull()

  ),
  
  # Tocilizumab arm
  
  "Events " = c(
    recovery_not_using %>% 
      select(events_toci) %>% pull(),
    recovery_using %>% 
      select(events_toci) %>% pull(),
    
   recovery_simple_oxygen %>% 
      select(events_toci) %>% pull(),
    recovery_noninvasive %>% 
      select(events_toci) %>% pull(),
   recovery_invasive %>% 
      select(events_toci) %>% pull()

  ),
  
  "Total " = c(
    
    recovery_not_using %>% 
      select(total_toci) %>% pull(),
    recovery_using %>% 
      select(total_toci) %>% pull(),
    
    recovery_simple_oxygen %>% 
      select(total_toci) %>% pull(),
    recovery_noninvasive %>% 
      select(total_toci) %>% pull(),
    recovery_invasive %>% 
      select(total_toci) %>% pull()

  ),

  "Risk (%)   " = c(
    
    recovery_not_using %>% 
      summarise(n = round(100*sum(events_toci)/sum(total_toci),1)) %>%
    pull(),
    recovery_using %>% 
      summarise(n = round(100*sum(events_toci)/sum(total_toci),1)) %>%
    pull(),
    
    recovery_simple_oxygen %>% 
      summarise(n = round(100*sum(events_toci)/sum(total_toci),1)) %>%
    pull(),
    recovery_noninvasive %>% 
      summarise(n = round(100*sum(events_toci)/sum(total_toci),1)) %>%
    pull(),
    recovery_invasive %>% 
      summarise(n = round(100*sum(events_toci)/sum(total_toci),1)) %>%
    pull()

  )
  
)

# Use kableExtra to have a nice table

# t1 = t %>% 
#   kbl(booktabs = T, align = 'c') %>% 
#   kable_styling(bootstrap_options = "striped", full_width = T) %>%
#   column_spec(c(5,8,11,14), bold = T) %>% 
#   add_header_above(c(" " = 2,
#                      "Control" = 3,
#                      "Tocilizumab" = 3,
#                      "Control" = 3,
#                      "Tocilizumab" = 3)) %>% 
#   add_header_above(c(" " = 1,
#                      "Priors" = 7,
#                      "RECOVERY" = 6)) %>% 
#   pack_rows("Use of corticosteroids", 1, 2) %>% 
#   pack_rows("Respiratory support", 3, 5)
# 
# t1

# t1 %>% 
#   save_kable(here("final_analyses", "output", "plots", "appendix", # File path
#                   "STable13_discharge_table1.pdf")) # File name

Prior and RECOVERY distributions

# Corticosteroids subgroups

# Create data frame

# Not using corticosteroids subgroup

distributions_not =
  posteriors_not_using_discharge %>%
  filter(type == "evidence-based") %>%
  summarise(
    RECOVERY_not = list(rnorm(10e4,
      mean = data.mean,
      sd = data.sd
    )),
    Prior_not = list(rnorm(10e4,
      mean = prior.mean,
      sd = prior.sd
    ))
  ) %>%
  unnest(RECOVERY_not:Prior_not)

distributions_not = bind_cols(
  distributions_not,
# Bind corresponding posterior distribution data
samples_posteriors_not_using_discharge %>% 
  filter(`Underlying_Prior` == "Evidence-based") %>% 
  select(`log(OR)`) %>% 
  rename("Posterior_not" = `log(OR)`)
)


# Using corticosteroids subgroup

distributions_u =
  posteriors_using_discharge %>%
  filter(type == "evidence-based") %>%
  summarise(
    RECOVERY_u = list(rnorm(10e4,
      mean = data.mean,
      sd = data.sd
    )),
    Prior_u = list(rnorm(10e4,
      mean = prior.mean,
      sd = prior.sd
    ))
  ) %>%
  unnest(RECOVERY_u:Prior_u)

distributions_u = bind_cols(
  distributions_u,
# Bind corresponding posterior distribution data
samples_posteriors_using_discharge %>% 
  filter(`Underlying_Prior` == "Evidence-based") %>% 
  select(`log(OR)`) %>% 
  rename("Posterior_u" = `log(OR)`)
)

# Bind altogether

distributions_cortico = bind_cols(
  distributions_not,
  distributions_u
  
)

arrow_labels = c("Tocilizumab Harm", "Tocilizumab Benefit")

xlim = c(0, 20)

xlab_df = data.frame(text = arrow_labels,
                          x = c(0.4,19),
                          y = c(0, 0),
                          hjust = c(0, 1))

a_small_amount = abs(xlim[1] - xlim[2])/20

null_line_at = 5

arrow_df = data.frame(id = c(1,2),
                      xstart = c(null_line_at - a_small_amount,
                                      null_line_at + a_small_amount),
                      xend = c(xlim[1] + a_small_amount, xlim[2] - a_small_amount),
                      y = c(1, 1))

arrows_plot = ggplot() +
      geom_segment(data = arrow_df,
                   aes(x = .data$xstart,
                       xend = .data$xend,
                       y = .data$y,
                       yend = .data$y),
                   arrow = arrow(angle = 15, type = "closed", length = grid::unit(0.1, "in"))) +
  geom_text(data = xlab_df,
            aes(x = .data$x,
                y = .data$y,
                label = .data$text,
                hjust = .data$hjust), size = 3.5) +
  scale_y_continuous(expand = c(0,0), limits = c(-0.5, 1.75)) +
  scale_x_continuous(expand = c(0,0), limits = xlim) +
  theme(panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.background = element_rect(fill = "transparent"),
        legend.box.background = element_rect(fill = "transparent"),
        panel.border = element_blank(),
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.line.y = element_blank(),
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.line.x = element_blank())



# Plot!

p1 = distributions_cortico %>%
  # make it tidy for ggplot
  pivot_longer(
    cols = c(RECOVERY_not:Posterior_u),
    names_to = "dist", values_to = "samples"
  ) %>%
  mutate(
    # Create new column to be able to use facet at the end
    subgroup = case_when(
      str_detect(dist, "_u") ~ "Using corticosteroids",
      str_detect(dist, "_not") ~ "Not using corticosteroids"
    ),

    # Set order of subgroups
    subgroup = factor(subgroup,
      levels = c(
        "Not using corticosteroids",
        "Using corticosteroids"
      )
    ),

    # Remove "_u" and "_not"
    dist = case_when(
      str_detect(dist, "RECOVERY") ~ "RECOVERY",
      str_detect(dist, "Prior") ~ "Prior",
      str_detect(dist, "Posterior") ~ "Posterior"
    ),
    # Set order of distributions
    dist = factor(dist,
      levels = c(
        "Posterior",
        "RECOVERY",
        "Prior"
      )
    )
  ) %>% 
  # exp() to convert log-odds ratio into odds ratio
  ggplot(aes(
    x = exp(samples), y = dist,
    fill = dist
  )) +

  # https://mjskay.github.io/ggdist/articles/slabinterval.html
  stat_halfeye(position = "dodge",
               point_interval = median_qi, # Quantile interval 
               .width = c(0.95)) +
  geom_vline(xintercept = 1, linetype = 2, size = 0.6) +
  scale_fill_manual(' ',
                    breaks=c("Prior","RECOVERY", "Posterior"),
                    values = c(
    "gray70", # Prior
    "#67B8D5", # RECOVERY
    "#F8D08D" # Posterior
  )) +
  labs(
    x = "\nOdds Ratio",
    y = " "
  ) +
  scale_x_continuous(breaks = seq(from = 0, to = 3, 0.5)) +
  coord_cartesian(xlim = c(0.25, 3),
                  clip = 'off') +
  scale_y_discrete(expand = c(0, 0.5)) +
  theme(
    strip.background = element_rect(fill = "#E4E6E7"),
    strip.text.x = element_text(size = 12),
    axis.ticks.x = element_blank(),
    axis.ticks.y = element_blank(),
    axis.text.x = element_text(size = 11),
    axis.text.y = element_blank(),
    panel.background = element_blank(),
    panel.grid.major.x = element_line(color = "gray80", size = 0.3),
    legend.position = 'right',
    legend.text = element_text(size=12),
    legend.key= element_blank(),
    plot.margin = margin(20, 0, 70, 25)
  ) +
  facet_wrap(~subgroup, ncol = 1)

## Respiratory support

# Create data frame

# Simple oxygen subgroup

distributions_s =
  posteriors_simple_discharge %>%
  filter(type == "evidence-based") %>%
  summarise(
    RECOVERY_s = list(rnorm(10e4,
      mean = data.mean,
      sd = data.sd
    )),
    Prior_s = list(rnorm(10e4,
      mean = prior.mean,
      sd = prior.sd
    ))
  ) %>%
  unnest(RECOVERY_s:Prior_s)

distributions_s = bind_cols(
  distributions_s,
# Bind corresponding posterior distribution data
samples_posteriors_simple_discharge %>% 
  filter(`Underlying_Prior` == "Evidence-based") %>% 
  select(`log(OR)`) %>% 
  rename("Posterior_s" = `log(OR)`)
)

# Non-invasive ventilation

distributions_n =
  posteriors_noninvasive_discharge %>%
  filter(type == "evidence-based") %>%
  summarise(
    RECOVERY_n = list(rnorm(10e4,
      mean = data.mean,
      sd = data.sd
    )),
    Prior_n = list(rnorm(10e4,
      mean = prior.mean,
      sd = prior.sd
    ))
  ) %>%
  unnest(RECOVERY_n:Prior_n)

distributions_n = bind_cols(
  distributions_n,
# Bind corresponding posterior distribution data
samples_posteriors_noninvasive_discharge %>% 
  filter(`Underlying_Prior` == "Evidence-based") %>% 
  select(`log(OR)`) %>% 
  rename("Posterior_n" = `log(OR)`)
)

# Invasive mechanical ventilation subgroup

distributions_i =
  posteriors_invasive_discharge %>%
  filter(type == "evidence-based") %>%
  summarise(
    RECOVERY_i = list(rnorm(10e4,
      mean = data.mean,
      sd = data.sd
    )),
    Prior_i = list(rnorm(10e4,
      mean = prior.mean,
      sd = prior.sd
    ))
  ) %>%
  unnest(RECOVERY_i:Prior_i)

distributions_i = bind_cols(
  distributions_i,
# Bind corresponding posterior distribution data
samples_posteriors_invasive_discharge %>% 
  filter(`Underlying_Prior` == "Evidence-based") %>% 
  select(`log(OR)`) %>% 
  rename("Posterior_i" = `log(OR)`)
)

# Bind altogether

distributions_respiratory = bind_cols(
  distributions_i,
  distributions_n,
  distributions_s
)

xlab_df = data.frame(text = arrow_labels,
                          x = c(0.4,19),
                          y = c(0, 0),
                          hjust = c(0, 1))

a_small_amount = abs(xlim[1] - xlim[2])/20

null_line_at = 5

arrow_df = data.frame(id = c(1,2),
                      xstart = c(null_line_at - a_small_amount,
                                      null_line_at + a_small_amount),
                      xend = c(xlim[1] + a_small_amount, xlim[2] - a_small_amount),
                      y = c(1, 1))

arrows_plot2 = ggplot() +
      geom_segment(data = arrow_df,
                   aes(x = .data$xstart,
                       xend = .data$xend,
                       y = .data$y,
                       yend = .data$y),
                   arrow = arrow(angle = 15, type = "closed", length = grid::unit(0.1, "in"))) +
  geom_text(data = xlab_df,
            aes(x = .data$x,
                y = .data$y,
                label = .data$text,
                hjust = .data$hjust), size = 3.5) +
  scale_y_continuous(expand = c(0,0), limits = c(-0.5, 1.75)) +
  scale_x_continuous(expand = c(0,0), limits = xlim) +
  theme(panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.background = element_rect(fill = "transparent"),
        legend.box.background = element_rect(fill = "transparent"),
        panel.border = element_blank(),
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.line.y = element_blank(),
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.line.x = element_blank())

# Plot!

p2 = distributions_respiratory %>%
  # make it tidy for ggplot
  pivot_longer(
    cols = c(RECOVERY_i:Prior_s),
    names_to = "dist", values_to = "samples"
  ) %>%
  mutate(
    # Create new column to be able to use facet at the end
    subgroup = case_when(
      str_detect(dist, "_s") ~ "Simple oxygen only",
      str_detect(dist, "_n") ~ "Non-invasive ventilation",
      str_detect(dist, "_i") ~ "Invasive mechanical ventilation"
    ),

    # Set order of subgroups
    subgroup = factor(subgroup,
      levels = c(
        "Simple oxygen only",
        "Non-invasive ventilation",
        "Invasive mechanical ventilation"
      )
    ),

    # Remove "_i", "_n" and "_s"
    dist = case_when(
      str_detect(dist, "RECOVERY") ~ "RECOVERY",
      str_detect(dist, "Prior") ~ "Prior",
      str_detect(dist, "Posterior") ~ "Posterior"
    ),
    # Set order of distributions
    dist = factor(dist,
      levels = c(
        "Posterior",
        "RECOVERY",
        "Prior"
      )
    )
  ) %>%
  # exp() to convert log-odds ratio into odds ratio
  ggplot(aes(
    x = exp(samples), y = dist,
    fill = dist
  )) +

  # https://mjskay.github.io/ggdist/articles/slabinterval.html
  stat_halfeye(position = "dodge",
               point_interval = median_qi, # Quantile interval
               .width = c(0.95)) +
  scale_fill_manual(values = c(
    "#F8D08D", # Posterior
    "#67B8D5", # RECOVERY
    "gray70" # Priors
    
  )) +
  geom_vline(xintercept = 1, linetype = 2, size = 0.6) +
  labs(
    x = "\nOdds Ratio",
    y = " "
  ) +
  scale_x_continuous(breaks = seq(0.5, 3.5, 0.5)) +
  coord_cartesian(xlim = c(0.8, 3.2)) +
  scale_y_discrete(expand = c(0, 0.5)) +
  theme(
    strip.background = element_rect(fill = "#E4E6E7"),
    strip.text.x = element_text(size = 12),
    axis.ticks.x = element_blank(),
    axis.ticks.y = element_blank(),
    axis.text.x = element_text(size = 11),
    axis.text.y = element_blank(),
    panel.background = element_blank(),
    panel.grid.major.x = element_line(color = "gray80", size = 0.3),
    legend.position = 'none',
    plot.margin = margin(20, 15, 70, 0)
  ) +
  facet_wrap(~subgroup, ncol = 1)

p1_overlay = p1 + inset_element(arrows_plot,
                                ignore_tag = TRUE,
                                align_to = "full",
                                left = unit(1.5, 'cm'),
                                bottom = unit(0.5, 'cm'),
                                right = unit(11.3, 'cm'),
                                top = unit(2.5, 'cm'))
# Plot everything

p1_overlay +
  p2 +
  inset_element(arrows_plot2,
                ignore_tag = TRUE,
                align_to = "full",
                left = unit(0, 'cm'),
                bottom = unit(0.5, 'cm'),
                right = unit(10.3, 'cm'),
                top = unit(2.5, 'cm')) +
  plot_annotation(tag_levels = 'A')

Point estimates depict the median and interval bars depict 95th quantiles.

# ggsave(width = 10,
#        height = 6,
#        here("final_analyses", "output", "plots", "appendix", # File path
#                  "SFigure13_discharge_figure1.pdf")) # File name

Here are the 95% quantile intervals in odds ratio for each distribution from the figure above.

t = distributions_cortico %>% 
  # make it tidy to be able to use group_by()
  pivot_longer(
    cols = c(RECOVERY_not:Posterior_u),
    names_to = "Subgroup", values_to = "samples"
  ) %>% 
  # Group by distribution so we can apply median_hdi
  group_by(Subgroup) %>% 
  # Calculate the median and 95% HDI in log-odds ratio
  ggdist::median_hdi(.width = 0.95) %>% 
  # Select only the relevant columns
  select(Subgroup:.upper) %>% 
  # Exponentiate to transform log-odds ratio into odds ratio
  mutate(across(samples:.upper, ~exp(.))) %>%
  mutate(across(samples:.upper, ~round(., 2))) %>%
  # Set order
  mutate(Subgroup = factor(Subgroup,
                           levels = c(
                             "Prior_not", "RECOVERY_not", "Posterior_not",
                             "Prior_u", "RECOVERY_u", "Posterior_u"
                           ))) %>% 
  arrange(Subgroup) %>%
  rename("Median" = samples,
         "Upper limit" = .upper,
         "Lower limit" = .lower)

# https://stackoverflow.com/questions/26548495/reorder-rows-using-custom-order
# create a vector with letters in the desired order
x = c("Prior_s", "RECOVERY_s", "Posterior_s",
      "Prior_n", "RECOVERY_n", "Posterior_n",
      "Prior_i", "RECOVERY_i", "Posterior_i")

t1 = distributions_respiratory %>% 
  # make it tidy to be able to use group_by()
  pivot_longer(
    cols = c(RECOVERY_i:Posterior_s),
    names_to = "Subgroup", values_to = "samples"
  ) %>% 
  # Group by distribution so we can apply median_hdi
  group_by(Subgroup) %>% 
  # Calculate the median and 95% HDI in log-odds ratio
  ggdist::median_hdi(.width = 0.95) %>% 
  # Select only the relevant columns
  select(Subgroup:.upper) %>% 
  # Exponentiate to transform log-odds ratio into odds ratio
  mutate(across(samples:.upper, ~exp(.))) %>% 
  mutate(across(samples:.upper, ~round(., 2))) %>%
  # Define order
  mutate(Subgroup = factor(Subgroup, levels = x)) %>% 
  arrange(Subgroup) %>%
  rename("Median" = samples,
         "Upper limit" = .upper,
         "Lower limit" = .lower)

tfinal = bind_rows(t, t1) %>% 
  mutate(
    Subgroup = case_when(
      str_detect(Subgroup, "Prior") ~ "Prior",
      str_detect(Subgroup, "RECOVERY") ~ "RECOVERY",
      str_detect(Subgroup, "Posterior") ~ "Posterior"
    )) %>% 
  rename("Distribution" = Subgroup) %>% 
  gt() %>%
  cols_width(
    Distribution ~ px(230)
  ) %>% 
  tab_row_group(
    label = "Invasive mechanical ventilation",
    rows = 13:15
  ) %>% 
  tab_row_group(
    label = "Non-invasive ventilation",
    rows = 10:12
  ) %>% 
  tab_row_group(
    label = "Simple oxygen only",
    rows = 7:9
  ) %>%
  tab_row_group(
    label = "Using corticosteroids",
    rows = 4:6
  ) %>% 
  tab_row_group(
    label = "Not using corticosteroids",
    rows = 1:3
  ) %>% 
  cols_align(
    align = "center",
    columns = everything()
  ) 

tfinal

Distribution	Median	Lower limit	Upper limit
Not using corticosteroids
Prior	1.66	0.94	2.96
RECOVERY	0.98	0.74	1.32
Posterior	1.10	0.84	1.42
Using corticosteroids
Prior	0.78	0.44	1.42
RECOVERY	1.42	1.24	1.62
Posterior	1.37	1.20	1.57
Simple oxygen only
Prior	1.59	0.83	2.94
RECOVERY	1.37	1.13	1.68
Posterior	1.39	1.15	1.68
Non-invasive ventilation
Prior	1.29	0.68	2.44
RECOVERY	1.34	1.10	1.62
Posterior	1.33	1.11	1.60
Invasive mechanical ventilation
Prior	1.45	0.70	3.03
RECOVERY	1.27	0.81	1.93
Posterior	1.31	0.90	1.92

# tfinal %>% 
#   gtsave(here("final_analyses", "output", "plots", "appendix", # File path
#                   "STable9_HDI_intervals_figureS13.png")) # File name

Posteriors distributions

# Create data frame

posterior_not_using = 
  samples_posteriors_not_using_discharge %>% 
  filter(`Underlying_Prior` == "Evidence-based") %>% 
  select(`RD`) %>% 
  rename("RD_not" = `RD`)

posterior_using = 
  samples_posteriors_using_discharge %>% 
  filter(`Underlying_Prior` == "Evidence-based") %>% 
  select(`RD`) %>% 
  rename("RD_u" = `RD`)

distributions_cortico_posterior =
  bind_cols(posterior_not_using,
            posterior_using)

# Prepare arrows under axis
# Adapted from https://github.com/rdboyes/forester/blob/master/R/forester.R


arrow_labels = c("Tocilizumab Benefit", "Tocilizumab Harm")

xlim = c(0, 20)

xlab_df = data.frame(text = arrow_labels,
                          x = c(2,19.5),
                          y = c(0, 0),
                          hjust = c(0, 1))

a_small_amount = abs(xlim[1] - xlim[2])/35

null_line_at = 14.5

arrow_df = data.frame(id = c(1,2),
                      xstart = c(null_line_at - a_small_amount,
                                      null_line_at + a_small_amount),
                      xend = c(xlim[1] + a_small_amount, xlim[2] - a_small_amount),
                      y = c(1, 1))

arrows_plot = ggplot() +
      geom_segment(data = arrow_df,
                   aes(x = .data$xstart,
                       xend = .data$xend,
                       y = .data$y,
                       yend = .data$y),
                   arrow = arrow(angle = 15, type = "closed", length = grid::unit(0.1, "in"))) +
  geom_text(data = xlab_df,
            aes(x = .data$x,
                y = .data$y,
                label = .data$text,
                hjust = .data$hjust), size = 3.5) +
  scale_y_continuous(expand = c(0,0), limits = c(-0.5, 1.75)) +
  scale_x_continuous(expand = c(0,0), limits = xlim) +
  theme(panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.background = element_rect(fill = "transparent"),
        legend.box.background = element_rect(fill = "transparent"),
        panel.border = element_blank(),
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.line.y = element_blank(),
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.line.x = element_blank())

p1 = distributions_cortico_posterior %>%
  # make it tidy for ggplot
  pivot_longer(
    cols = c(RD_not:RD_u),
    names_to = "dist", values_to = "samples"
  ) %>%
  mutate(
    # Create new column to be able to use facet at the end
    subgroup = case_when(
      str_detect(dist, "_u") ~ "Using corticosteroids",
      str_detect(dist, "_not") ~ "Not using corticosteroids"
    ),

    # Set order of subgroups
    subgroup = factor(subgroup,
      levels = c(
        "Not using corticosteroids",
        "Using corticosteroids"
      )
    )) %>%

  # https://mjskay.github.io/ggdist/articles/slabinterval.html
  ggplot(aes(
    # Multiply by 100 to report in percentage
    x = 100 * samples, y = dist,
    fill = subgroup,
   # https://github.com/mjskay/ggdist/issues/71
    fill_ramp = stat(x < 0))
  ) +
  # https://mjskay.github.io/ggdist/articles/slabinterval.html
  stat_halfeye(point_interval = median_hdi, # Highest density interval
               .width = c(0.8, 0.95),
               # https://github.com/mjskay/ggdist/issues/70
               interval_size_range = c(1.1,1.9)) +
  scale_fill_ramp_discrete(from = "gray85", range = c(0,1)) +
  scale_fill_manual(values = c("#9EA45B", "#A65041")) +
  geom_vline(xintercept = 0, linetype = 2, size = 0.6) +
  labs(
    x = "\nRisk Difference (%)",
    y = " "
  ) +
  scale_x_continuous(breaks = seq(from = -15, to = 5, 2.5)) +
  coord_cartesian(xlim = c(-15, 5)) +
  scale_y_discrete(expand = c(0, 0.5)) +
  theme(
    strip.background = element_rect(fill = "#E4E6E7"),
    strip.text.x = element_text(size = 12),
    axis.ticks.x = element_blank(),
    axis.ticks.y = element_blank(),
    axis.text.x = element_text(size = 11),
    axis.text.y = element_blank(),
    panel.background = element_blank(),
    panel.grid.major.x = element_line(color = "gray80", size = 0.3),
    legend.position = 'none',
    plot.margin = margin(10, 10, 50, 10)
  ) +
  facet_wrap(~subgroup, ncol = 1, scales = "free_y")

# Create data frame

# Simple oxygen only

posterior_s = 
  samples_posteriors_simple_discharge %>% 
  filter(`Underlying_Prior` == "Evidence-based") %>% 
  select(`RD`) %>% 
  rename("RD_s" = `RD`)

# Non-invasive ventilation

posterior_n = 
  samples_posteriors_noninvasive_discharge %>% 
  filter(`Underlying_Prior` == "Evidence-based") %>% 
  select(`RD`) %>% 
  rename("RD_n" = `RD`)

# Invasive mechanical ventilation

posterior_i = 
  samples_posteriors_invasive_discharge %>% 
  filter(`Underlying_Prior` == "Evidence-based") %>% 
  select(`RD`) %>% 
  rename("RD_i" = `RD`)

# Bind altogether

distributions_respiratory_posterior = bind_cols(
  posterior_i,
  posterior_n,
  posterior_s
)

p3 = distributions_respiratory_posterior %>%
  # make it tidy for ggplot
  pivot_longer(
    cols = c(RD_i:RD_s),
    names_to = "dist", values_to = "samples"
  ) %>%
  mutate(
    # Create new column to be able to use facet at the end
    subgroup = case_when(
      str_detect(dist, "_s") ~ "Simple oxygen only",
      str_detect(dist, "_n") ~ "Non-invasive ventilation",
      str_detect(dist, "_i") ~ "Invasive mechanical ventilation"
    ),

    # Set order of subgroups
    subgroup = factor(subgroup,
      levels = c(
        "Simple oxygen only",
        "Non-invasive ventilation",
        "Invasive mechanical ventilation"
      ))
    ) %>%
  # Plot!
  ggplot(aes(
    # Multiply by 100 to report in percentage
    x = 100 * samples, y = dist,
    fill = subgroup,
    # https://github.com/mjskay/ggdist/issues/71
    fill_ramp = stat(x < 0))
  ) +

  # https://mjskay.github.io/ggdist/articles/slabinterval.html
  stat_halfeye(point_interval = median_hdi, # Highest density interval
               .width = c(0.8, 0.95),
               # https://github.com/mjskay/ggdist/issues/70
               interval_size_range = c(1.1,1.9)) +
  scale_fill_ramp_discrete(from = "gray85", range = c(0,1)) +
  scale_fill_manual(values = c("#9D75B1", # Simple oxygen
                               "#CE9F5B", # Non-invasive
                               "#5891BF") # Invasive
                    ) +
  geom_vline(xintercept = 0, linetype = 2, size = 0.6) +
  labs(
    x = "\nRisk Difference (%)",
    y = " "
  ) +
  scale_x_continuous(breaks = seq(from = -15, to = 5, 2.5)) +
  coord_cartesian(xlim = c(-15, 5)) +
  scale_y_discrete(expand = c(0, 0.5)) +
  theme(
    strip.background = element_rect(fill = "#E4E6E7"),
    strip.text.x = element_text(size = 12),
    axis.ticks.x = element_blank(),
    axis.ticks.y = element_blank(),
    axis.text.x = element_text(size = 11),
    axis.text.y = element_blank(),
    panel.background = element_blank(),
    panel.grid.major.x = element_line(color = "gray80", size = 0.3),
    legend.position = 'none',
    plot.margin = margin(20, 10, 50, 10)
  ) +
  facet_wrap(~ subgroup, ncol = 1, scales = "free_y")

p4 = distributions_respiratory_posterior %>%
  # make it tidy for ggplot
  pivot_longer(
    cols = c(RD_i:RD_s),
    names_to = "dist", values_to = "samples"
  ) %>%
  mutate(
    # Create new column to be able to use facet at the end
    subgroup = case_when(
      str_detect(dist, "_s") ~ "Simple oxygen only",
      str_detect(dist, "_n") ~ "Non-invasive ventilation",
      str_detect(dist, "_i") ~ "Invasive mechanical ventilation"
    ),

    # Set order of subgroups
    subgroup = factor(subgroup,
      levels = c(
        "Simple oxygen only",
        "Non-invasive ventilation",
        "Invasive mechanical ventilation"
      ))
    ) %>%
# Multiply by 100 to report in percentage
ggplot(aes(100*samples, colour = subgroup)) +
  stat_ecdf(geom = "step", size = 1.2) +
  scale_color_manual(' ',
                    values = c("#9D75B1", # Simple oxygen
                               "#CE9F5B", # Non-invasive
                               "#5891BF") # Invasive
                    ) +
  geom_vline(xintercept = 0, linetype = 2, size = 0.6) +
    labs(
      x = "\nRisk Difference (%)",
      y = "Probability RD \u2264 X\n"
      ) +
    scale_x_continuous(
      breaks = seq(from = -15, to = 10, 1),
      limits = c(-15, 10)) +
    scale_y_continuous(
      breaks = seq(from = 0, to = 1, .1),
      expand = c(0, .03)
    ) +
  theme(
      axis.ticks.x = element_blank(),
      panel.background = element_blank(),
      panel.grid.major.x = element_line(color = "gray80", size = 0.3),
      panel.grid.major.y = element_line(color = "gray80", size = 0.3),
      legend.position = 'none',
      plot.margin = margin(10, 10, 50, 10)
    ) +
  coord_cartesian(x = c(-12.5,0.1))

p1_overlay = p1 + inset_element(arrows_plot,
                                ignore_tag = TRUE,
                                align_to = "full",
                                left = unit(0.65, 'cm'),
                                bottom = unit(0, 'cm'),
                                right = unit(12.3, 'cm'),
                                top = unit(1.9, 'cm'))
# Plot everything



(p1_overlay + p2) / (p3 + 
                      inset_element(arrows_plot,
                                ignore_tag = TRUE,
                                align_to = "full",
                                left = unit(0.65, 'cm'),
                                bottom = unit(0, 'cm'),
                                right = unit(12.3, 'cm'),
                                top = unit(1.9, 'cm')) + 
                       p4) +
  plot_annotation(tag_levels = 'A')

Panels A and C: Point estimates depict the median. Interval bars depict the 80 and 95% highest density intervals. Panels B and D: Cumulative posterior distributions correspond to the probabilities that the risk difference (RD) is less than or equal to the effect size on the X‐axis.

# ggsave(width = 10,
#        height = 8.5,
#  #https://stackoverflow.com/questions/12768176/unicode-characters-in-ggplot2-pdf-output
#        device=cairo_pdf, 
#        here("final_analyses", "output", "plots", "appendix", # File path
#                  "SFigure14_discharge_figure2.pdf")) # File name

multiply100 = function(x){x*100}

t1 = distributions_cortico_posterior %>%
  # make it tidy to be able to use group_by()
  pivot_longer(
    cols = c(RD_not:RD_u),
    names_to = "Subgroup", values_to = "samples"
  ) %>% 
  group_by(Subgroup) %>% 
  summarise("Pr(\u2264 0%)" = (sum( samples <= 0))/ n(),
            "Pr(\u2264 -1%)" = (sum( samples <= -0.01))/ n(),
            "Pr(\u2264 -2%)" = (sum( samples <= -0.02))/ n(),
            "Pr(\u2264 -5%)" = (sum( samples <= -0.05))/ n()) %>% 
  mutate(across(-Subgroup, ~multiply100(.))) %>% 
  mutate(across(-Subgroup, ~round(.,1))) %>% 
  mutate(
    Subgroup = case_when(
      str_detect(Subgroup, "_u") ~ "Using",
      str_detect(Subgroup, "_not") ~ "Not using"
    ))

multiply100 = function(x){x*100}
t2 = distributions_cortico_posterior %>% 
  # make it tidy to be able to use group_by()
  pivot_longer(
    cols = c(RD_not:RD_u),
    names_to = "Subgroup", values_to = "samples"
  ) %>% 
  # Group by distribution so we can apply median_hdi
  group_by(Subgroup) %>% 
  # Calculate the median and 95% quantile interval in log-odds ratio
  ggdist::median_hdi(.width = 0.95) %>% 
  # Select only the relevant columns
  select(Subgroup:.upper) %>%
  mutate(across(samples:.upper, ~multiply100(.))) %>%
  mutate(across(samples:.upper, ~round(., 2))) %>%
  rename("Median" = samples,
         "Upper limit" = .upper,
         "Lower limit" = .lower) %>% 
  mutate(
    Subgroup = case_when(
      str_detect(Subgroup, "_u") ~ "Using",
      str_detect(Subgroup, "_not") ~ "Not using"
    ))

multiply100 = function(x){x*100}

t3 = distributions_respiratory_posterior %>%
  # make it tidy to be able to use group_by()
  pivot_longer(
    cols = c(RD_i:RD_s),
    names_to = "Subgroup", values_to = "samples"
  ) %>% 
  group_by(Subgroup) %>% 
  summarise("Pr(\u2264 0%)" = (sum( samples <= 0))/ n(),
            "Pr(\u2264 -1%)" = (sum( samples <= -0.01))/ n(),
            "Pr(\u2264 -2%)" = (sum( samples <= -0.02))/ n(),
            "Pr(\u2264 -5%)" = (sum( samples <= -0.05))/ n()) %>% 
  mutate(across(-Subgroup, ~multiply100(.))) %>% 
  mutate(across(-Subgroup, ~round(.,1))) %>%
  arrange(desc(Subgroup)) %>% 
  mutate(
    Subgroup = case_when(
      str_detect(Subgroup, "_s") ~ "Simple oxygen only",
      str_detect(Subgroup, "_n") ~ "Non-invasive ventilation",
      str_detect(Subgroup, "_i") ~ "Invasive mechanical ventilation"
    ))

multiply100 = function(x){x*100}

t4 = distributions_respiratory_posterior %>% 
  # make it tidy to be able to use group_by()
  pivot_longer(
    cols = c(RD_i:RD_s),
    names_to = "Subgroup", values_to = "samples"
  ) %>% 
  # Group by distribution so we can apply median_hdi
  group_by(Subgroup) %>% 
  # Calculate the median and 95% quantile interval in log-odds ratio
  ggdist::median_hdi(.width = 0.95) %>% 
  # Select only the relevant columns
  select(Subgroup:.upper) %>%
  mutate(across(samples:.upper, ~multiply100(.))) %>%
  mutate(across(samples:.upper, ~round(., 2))) %>%
  rename("Median" = samples,
         "Upper limit" = .upper,
         "Lower limit" = .lower) %>% 
  arrange(desc(Subgroup)) %>% 
  mutate(
    Subgroup = case_when(
      str_detect(Subgroup, "_s") ~ "Simple oxygen only",
      str_detect(Subgroup, "_n") ~ "Non-invasive ventilation",
      str_detect(Subgroup, "_i") ~ "Invasive mechanical ventilation"
    ))

Here are the posterior probabilities for each posterior distribution.

bind_rows(t1,t3) %>% 
  kbl(booktabs = T, align = 'c') %>% 
  kable_styling(bootstrap_options = "striped", full_width = F) %>%
  pack_rows("Use of corticosteroids", 1, 2) %>% 
  pack_rows("Respiratory support", 3, 5) %>% 
  add_header_above(c(" " = 1,
                     "Probability of Risk Difference \u2265 X%" = 4))

	Probability of Risk Difference ≥ X%
Subgroup	Pr(≤ 0%)	Pr(≤ -1%)	Pr(≤ -2%)	Pr(≤ -5%)
Use of corticosteroids
Not using	75.6	65.2	53.6	20.8
Using	100.0	100.0	100.0	95.9
Respiratory support
Simple oxygen only	100.0	99.8	99.3	82.2
Non-invasive ventilation	99.9	99.6	98.6	81.8
Invasive mechanical ventilation	92.1	84.8	74.8	36.8

Here are the 95% HDI in risk difference of each posterior distribution.

bind_rows(t2,t4) %>% 
  kbl(booktabs = T, align = 'c') %>% 
  kable_styling(bootstrap_options = "striped", full_width = F) %>%
  pack_rows("Use of corticosteroids", 1, 2) %>% 
  pack_rows("Respiratory support", 3, 5) %>% 
  add_header_above(c(" " = 1,
                     "95% Highest Density Interval" = 3))

	95% Highest Density Interval
Subgroup	Median	Lower limit	Upper limit
Use of corticosteroids
Not using	-2.30	-8.77	4.16
Using	-7.87	-11.04	-4.63
Respiratory support
Simple oxygen only	-6.75	-10.27	-3.08
Non-invasive ventilation	-7.13	-11.74	-2.56
Invasive mechanical ventilation	-3.96	-10.33	1.58

Sensitivity analyses: different priors

Not using corticosteroids

Prior and likelihood distributions

distributions =
  tibble(
    
    "RECOVERY (likelihood)" = list(
      rnorm(
        10e4,
        mean = (posteriors_not_using_discharge %>%
                # Data.mean and data.sd are all the same in every row
                # So it doesn't matter what row I slice
                slice(1) %>% 
                pull(data.mean)),
        sd = (posteriors_not_using_discharge %>%
              slice(1) %>% 
              pull(data.sd))
    )),
    
    "Non-informative prior" = list(
      rnorm(
        10e4,
        mean = (posteriors_not_using_discharge %>%
                filter(type == "non-informative") %>%
                pull(prior.mean)),
        sd = (posteriors_not_using_discharge %>%
              filter(type == "non-informative") %>%
              pull(prior.sd))
    )),
    
    "Evidence-based prior" = list(
      rnorm(
        10e4,
        mean = (posteriors_not_using_discharge %>%
                filter(type == "evidence-based") %>%
                pull(prior.mean)),
        sd = (posteriors_not_using_discharge %>%
              filter(type == "evidence-based") %>%
              pull(prior.sd))
    )),
    
    "Skeptical prior" = list(
      rnorm(
        10e4,
        mean = (posteriors_not_using_discharge %>%
                filter(type == "skeptical") %>%
                pull(prior.mean)),
        sd = (posteriors_not_using_discharge %>%
              filter(type == "skeptical") %>%
              pull(prior.sd))
    )),
    
    "Optimistic prior" = list(
      rnorm(
        10e4,
        mean = (posteriors_not_using_discharge %>%
                filter(type == "optimistic") %>%
                pull(prior.mean)),
        sd = (posteriors_not_using_discharge %>%
              filter(type == "optimistic") %>%
              pull(prior.sd))
    )),
    
    "Pessimistic prior" = list(
      rnorm(
        10e4,
        mean = (posteriors_not_using_discharge %>%
                filter(type == "pessimistic") %>%
                pull(prior.mean)),
        sd = (posteriors_not_using_discharge %>%
              filter(type == "pessimistic") %>%
              pull(prior.sd))
    ))
  ) %>% 
  unnest(everything())

# Plot!

distributions %>%
  # make it tidy for ggplot
  pivot_longer(
    cols = c(everything()),
    names_to = "dist", values_to = "samples"
  ) %>%
  # Set order of priors
  mutate(dist = factor(dist, levels =
                         c("RECOVERY (likelihood)",
                           "Non-informative prior",
                           "Evidence-based prior",
                           "Skeptical prior",
                           "Optimistic prior",
                           "Pessimistic prior")
                       )
         ) %>%
  # exp() to transform log(OR) to OR
  ggplot(aes(
    x = exp(samples), y = dist,
    fill = dist
  )) +

  # https://mjskay.github.io/ggdist/articles/slabinterval.html
  stat_halfeye(position = "dodge",
               point_interval = median_qi, # Quantile interval 
               .width = c(0.95)) +
  geom_vline(xintercept = 1, linetype = 2, size = 0.6) +
  scale_fill_npg()+
  labs(
    x = "\nOdds Ratio",
    y = " "
  ) +
  scale_x_continuous(breaks = seq(from = 0, to = 3, 0.5)) +
  scale_y_discrete(expand = c(0, 0.5)) +
  theme(
    strip.background = element_rect(fill = "#E4E6E7"),
    strip.text.x = element_text(size = 12),
    axis.ticks.x = element_blank(),
    axis.ticks.y = element_blank(),
    axis.text.x = element_text(size = 11),
    axis.text.y = element_blank(),
    panel.background = element_blank(),
    panel.grid.major.x = element_line(color = "gray80", size = 0.3),
    legend.position = 'none',
    legend.text = element_text(size=12),
    legend.key= element_blank(),
    plot.margin = margin(20, 20, 20, 20)
  ) +
  coord_cartesian(x = c(0, 3)) +
  facet_wrap(~dist, ncol = 1, scales = "free_y")

Posterior distributions

# Prepare arrows under axis
# Adapted from https://github.com/rdboyes/forester/blob/master/R/forester.R


arrow_labels = c("Tocilizumab Benefit", "Tocilizumab Harm")

xlim = c(0, 20)

xlab_df = data.frame(text = arrow_labels,
                          x = c(2,18),
                          y = c(0, 0),
                          hjust = c(0, 1))

a_small_amount = abs(xlim[1] - xlim[2])/35

null_line_at = 10

arrow_df = data.frame(id = c(1,2),
                      xstart = c(null_line_at - a_small_amount,
                                      null_line_at + a_small_amount),
                      xend = c(xlim[1] + a_small_amount, xlim[2] - a_small_amount),
                      y = c(1, 1))

arrows_plot = ggplot() +
      geom_segment(data = arrow_df,
                   aes(x = .data$xstart,
                       xend = .data$xend,
                       y = .data$y,
                       yend = .data$y),
                   arrow = arrow(angle = 15, type = "closed", length = grid::unit(0.1, "in"))) +
  geom_text(data = xlab_df,
            aes(x = .data$x,
                y = .data$y,
                label = .data$text,
                hjust = .data$hjust), size = 3.5) +
  scale_y_continuous(expand = c(0,0), limits = c(-0.5, 1.75)) +
  scale_x_continuous(expand = c(0,0), limits = xlim) +
  theme(panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.background = element_rect(fill = "transparent"),
        legend.box.background = element_rect(fill = "transparent"),
        panel.border = element_blank(),
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.line.y = element_blank(),
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.line.x = element_blank())

level_order = c("Pessimistic",
                "Optimistic",
                "Skeptical",
                "Evidence-based",
                "Non-informative")

p1 = samples_posteriors_not_using_discharge %>%
  
  ggplot(aes(
    # Multiply by 100 to report in percentage
    x = 100 * RD,
    y = factor(`Underlying_Prior`, level = level_order),
    fill = `Underlying_Prior`)) +
  geom_vline(xintercept = seq(-12.5, 12.5, 2.5), linetype = 1, size = 0.3,
             color = "gray80") +

  # https://mjskay.github.io/ggdist/articles/slabinterval.html
    stat_halfeye(aes(fill = stat(x < 0)),
                 position = "dodge",
               point_interval = median_hdi, # Highest density interval
               .width = c(0.8, 0.95),
               # https://github.com/mjskay/ggdist/issues/70
               interval_size_range = c(1.1,1.9)
  ) +
  scale_fill_manual(values = c("gray85", "#9EA45B")) +
  labs(x = "\nRisk Difference (%)",
       y = "Underlying Prior\n"
  ) +
  scale_x_continuous(breaks = seq(from = -15, to = 15, 5)) +
  scale_y_discrete(expand = c(0, 0.3)) +
  theme(axis.title.x = element_text(size = 12),
    axis.title.y = element_text(size = 13.5),
    axis.text.y = element_text(size = 12),
    axis.text.x = element_text(size = 12),
    axis.ticks.x = element_blank(),
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    panel.grid.major.x = element_line(color = "gray80", size = 0.3),
    legend.position = "none",
    plot.title.position = 'plot',
    plot.margin = margin(20, 25, 40, 20)
  ) +
  coord_cartesian(xlim = c(-15, 15))

p1_overlay = p1 + inset_element(arrows_plot,
                                ignore_tag = TRUE,
                                align_to = "full",
                                left = unit(4.8, 'cm'),
                                bottom = unit(0, 'cm'),
                                right = unit(16.8, 'cm'),
                                top = unit(1.6, 'cm'))

p1_overlay

Point estimates depict the median. Interval bars depict the 80 and 95% highest density intervals.

# ggsave(width = 7,
#        height = 4.5,
#         here("final_analyses", "output", "plots", "appendix", # File path
#                   "SFigure15_discharge_posteriors_not_using.pdf")) # File name

t = tibble(
  
  "Underlying Prior" = c(
    "Non-informative",
    "Evidence-based",
    "Skeptical",
    "Optimistic",
    "Pessimistic"
  ),
  
  
  "Pr(< -5%)" = c(round((samples_posteriors_not_using_discharge %>%
                           filter(`Underlying_Prior` == "Non-informative") %>%
                           summarise(n = sum(RD < -0.05) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_not_using_discharge %>%
                           filter(`Underlying_Prior` == "Evidence-based") %>%
                           summarise(n = sum(RD < -0.05) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_not_using_discharge %>%
                           filter(`Underlying_Prior` == "Skeptical") %>%
                           summarise(n = sum(RD < -0.05) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_not_using_discharge %>%
                           filter(`Underlying_Prior` == "Optimistic") %>%
                           summarise(n = sum(RD < -0.05) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_not_using_discharge %>%
                           filter(`Underlying_Prior` == "Pessimistic") %>%
                           summarise(n = sum(RD < -0.05) / n()) %>%
                           pull()), 2) * 100
  ),
  
  "Pr(< -2%)" = c(round((samples_posteriors_not_using_discharge %>%
                           filter(`Underlying_Prior` == "Non-informative") %>%
                           summarise(n = sum(RD < -0.02) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_not_using_discharge %>%
                           filter(`Underlying_Prior` == "Evidence-based") %>%
                           summarise(n = sum(RD < -0.02) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_not_using_discharge %>%
                           filter(`Underlying_Prior` == "Skeptical") %>%
                           summarise(n = sum(RD < -0.02) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_not_using_discharge %>%
                           filter(`Underlying_Prior` == "Optimistic") %>%
                           summarise(n = sum(RD < -0.02) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_not_using_discharge %>%
                           filter(`Underlying_Prior` == "Pessimistic") %>%
                           summarise(n = sum(RD < -0.02) / n()) %>%
                           pull()), 2) * 100
  ),
  
  "Pr(< -1%)" = c(round((samples_posteriors_not_using_discharge %>%
                           filter(`Underlying_Prior` == "Non-informative") %>%
                           summarise(n = sum(RD < -0.01) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_not_using_discharge %>%
                           filter(`Underlying_Prior` == "Evidence-based") %>%
                           summarise(n = sum(RD < -0.01) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_not_using_discharge %>%
                           filter(`Underlying_Prior` == "Skeptical") %>%
                           summarise(n = sum(RD < -0.01) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_not_using_discharge %>%
                           filter(`Underlying_Prior` == "Optimistic") %>%
                           summarise(n = sum(RD < -0.01) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_not_using_discharge %>%
                           filter(`Underlying_Prior` == "Pessimistic") %>%
                           summarise(n = sum(RD < -0.01) / n()) %>%
                           pull()), 2) * 100
  ),
  
  
  "Pr(< 0%)" = c(round((samples_posteriors_not_using_discharge %>%
                          filter(`Underlying_Prior` == "Non-informative") %>%
                          summarise(n = sum(RD < 0) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_not_using_discharge %>%
                          filter(`Underlying_Prior` == "Evidence-based") %>%
                          summarise(n = sum(RD < 0) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_not_using_discharge %>%
                          filter(`Underlying_Prior` == "Skeptical") %>%
                          summarise(n = sum(RD < 0) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_not_using_discharge %>%
                          filter(`Underlying_Prior` == "Optimistic") %>%
                          summarise(n = sum(RD < 0) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_not_using_discharge %>%
                          filter(`Underlying_Prior` == "Pessimistic") %>%
                          summarise(n = sum(RD < 0) / n()) %>%
                          pull()), 2) * 100
  ),
  
  "Pr(> 1%)" = c(round((samples_posteriors_not_using_discharge %>%
                          filter(`Underlying_Prior` == "Non-informative") %>%
                          summarise(n = sum(RD > 0.01) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_not_using_discharge %>%
                          filter(`Underlying_Prior` == "Evidence-based") %>%
                          summarise(n = sum(RD > 0.01) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_not_using_discharge %>%
                          filter(`Underlying_Prior` == "Skeptical") %>%
                          summarise(n = sum(RD > 0.01) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_not_using_discharge %>%
                          filter(`Underlying_Prior` == "Optimistic") %>%
                          summarise(n = sum(RD > 0.01) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_not_using_discharge %>%
                          filter(`Underlying_Prior` == "Pessimistic") %>%
                          summarise(n = sum(RD > 0.01) / n()) %>%
                          pull()), 2) * 100
  ),
  
  "Pr(> 2%)" = c(round((samples_posteriors_not_using_discharge %>%
                          filter(`Underlying_Prior` == "Non-informative") %>%
                          summarise(n = sum(RD > 0.02) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_not_using_discharge %>%
                          filter(`Underlying_Prior` == "Evidence-based") %>%
                          summarise(n = sum(RD > 0.02) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_not_using_discharge %>%
                          filter(`Underlying_Prior` == "Skeptical") %>%
                          summarise(n = sum(RD > 0.02) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_not_using_discharge %>%
                          filter(`Underlying_Prior` == "Optimistic") %>%
                          summarise(n = sum(RD > 0.02) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_not_using_discharge %>%
                          filter(`Underlying_Prior` == "Pessimistic") %>%
                          summarise(n = sum(RD > 0.02) / n()) %>%
                          pull()), 2) * 100
  ),
  
  "Pr(> 5%)" = c(round((samples_posteriors_not_using_discharge %>%
                          filter(`Underlying_Prior` == "Non-informative") %>%
                          summarise(n = sum(RD > 0.05) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_not_using_discharge %>%
                          filter(`Underlying_Prior` == "Evidence-based") %>%
                          summarise(n = sum(RD > 0.05) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_not_using_discharge %>%
                          filter(`Underlying_Prior` == "Skeptical") %>%
                          summarise(n = sum(RD > 0.05) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_not_using_discharge %>%
                          filter(`Underlying_Prior` == "Optimistic") %>%
                          summarise(n = sum(RD > 0.05) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_not_using_discharge %>%
                          filter(`Underlying_Prior` == "Pessimistic") %>%
                          summarise(n = sum(RD > 0.05) / n()) %>%
                          pull()), 2) * 100
  )
) 

# https://stackoverflow.com/questions/26548495/reorder-rows-using-custom-order

level_order = c("Evidence-based",
                "Non-informative",
                "Skeptical",
                "Optimistic",
                "Pessimistic"
                )

t = t %>% 
  slice(match(level_order, `Underlying Prior`))

HDI for each posterior distribution

Here are the 95% HDI in risk difference of each distribution.

multiply100 = function(x){x*100}

t1 = samples_posteriors_not_using_discharge %>%
  select("Underlying_Prior", "RD") %>% 
  # Group by distribution so we can apply median_hdi
  group_by(`Underlying_Prior`) %>% 
  # Calculate the median and 95% quantile interval in log-odds ratio
  ggdist::median_hdi(.width = 0.95) %>% 
  # Select only the relevant columns
  # Arrange order
  arrange(factor(`Underlying_Prior`, levels = c(
    "Evidence-based",
    "Non-informative",
    "Skeptical",
    "Optimistic",
    "Pessimistic"
    ))) %>% 
  select("Underlying_Prior":.upper) %>%
  mutate(across(RD:.upper, ~multiply100(.))) %>%
  mutate(across(RD:.upper, ~round(., 2))) %>%
  rename("Underlying Prior" = Underlying_Prior,
         "Median" = RD,
         "Lower limit" = .lower,
         "Upper limit" = .upper) 

tfinal = left_join(t1, t)

## Joining, by = "Underlying Prior"

tfinal = tfinal %>% 
  gt() %>% 
  tab_spanner(label = "95% HDI",
              columns = c(2:4)) %>% 
  tab_spanner(label = "Propability of Benefit",
              columns = c("Pr(< -5%)", "Pr(< -2%)", "Pr(< -1%)", "Pr(< 0%)")) %>% 
  tab_spanner(label = "Probability of Harm",
              columns = c("Pr(> 1%)", "Pr(> 2%)", "Pr(> 5%)")) %>%
  cols_align(
    align = "center",
    columns = everything()
  ) 
  

tfinal

Underlying Prior	95% HDI			Propability of Benefit				Probability of Harm
	Median	Lower limit	Upper limit	Pr(< -5%)	Pr(< -2%)	Pr(< -1%)	Pr(< 0%)	Pr(> 1%)	Pr(> 2%)	Pr(> 5%)
Evidence-based	-2.30	-8.77	4.16	21	54	65	76	16	10	1
Non-informative	0.44	-6.90	7.51	7	26	35	45	44	34	11
Skeptical	0.33	-6.15	6.74	6	24	35	46	42	31	8
Optimistic	-0.82	-7.30	5.60	10	36	48	60	29	19	4
Pessimistic	1.46	-4.87	7.91	3	15	23	33	55	43	14

# tfinal %>% 
#   gtsave(here("final_analyses", "output", "plots", "appendix", # File path
#                   "STable10_discharge_hdi_probabilities_not_using.png"))

Using corticosteroids

Prior and likelihood distributions

distributions =
  tibble(
    
    "RECOVERY (likelihood)" = list(
      rnorm(
        10e4,
        mean = (posteriors_using_discharge %>%
                  # Data.mean and data.sd are all the same in every row
                  # So it doesn't matter what row I slice
                  slice(1) %>% 
                  pull(data.mean)),
        sd = (posteriors_using_discharge %>%
                slice(1) %>% 
                pull(data.sd))
      )),
    
    "Non-informative prior" = list(
      rnorm(
        10e4,
        mean = (posteriors_using_discharge %>%
                  filter(type == "non-informative") %>%
                  pull(prior.mean)),
        sd = (posteriors_using_discharge %>%
                filter(type == "non-informative") %>%
                pull(prior.sd))
      )),
    
    "Evidence-based prior" = list(
      rnorm(
        10e4,
        mean = (posteriors_using_discharge %>%
                  filter(type == "evidence-based") %>%
                  pull(prior.mean)),
        sd = (posteriors_using_discharge %>%
                filter(type == "evidence-based") %>%
                pull(prior.sd))
      )),
    
    "Skeptical prior" = list(
      rnorm(
        10e4,
        mean = (posteriors_using_discharge %>%
                  filter(type == "skeptical") %>%
                  pull(prior.mean)),
        sd = (posteriors_using_discharge %>%
                filter(type == "skeptical") %>%
                pull(prior.sd))
      )),
    
    "Optimistic prior" = list(
      rnorm(
        10e4,
        mean = (posteriors_using_discharge %>%
                  filter(type == "optimistic") %>%
                  pull(prior.mean)),
        sd = (posteriors_using_discharge %>%
                filter(type == "optimistic") %>%
                pull(prior.sd))
      )),
    
    "Pessimistic prior" = list(
      rnorm(
        10e4,
        mean = (posteriors_using_discharge %>%
                  filter(type == "pessimistic") %>%
                  pull(prior.mean)),
        sd = (posteriors_using_discharge %>%
                filter(type == "pessimistic") %>%
                pull(prior.sd))
      ))
  ) %>% 
  unnest(everything())

# Plot!

distributions %>%
  # make it tidy for ggplot
  pivot_longer(
    cols = c(everything()),
    names_to = "dist", values_to = "samples"
  ) %>%
  # Set order of priors
  mutate(dist = factor(dist, levels =
                         c("RECOVERY (likelihood)",
                           "Non-informative prior",
                           "Evidence-based prior",
                           "Skeptical prior",
                           "Optimistic prior",
                           "Pessimistic prior")
  )
  ) %>%
  # exp() to transform log(OR) to OR
  ggplot(aes(
    x = exp(samples), y = dist,
    fill = dist
  )) +
  
  # https://mjskay.github.io/ggdist/articles/slabinterval.html
  stat_halfeye(position = "dodge",
               point_interval = median_qi, # Quantile interval 
               .width = c(0.95)) +
  geom_vline(xintercept = 1, linetype = 2, size = 0.6) +
  scale_fill_npg()+
  labs(
    x = "\nOdds Ratio",
    y = " "
  ) +
  scale_x_continuous(breaks = seq(from = 0, to = 3, 0.5)) +
  scale_y_discrete(expand = c(0, 0.5)) +
  theme(
    strip.background = element_rect(fill = "#E4E6E7"),
    strip.text.x = element_text(size = 12),
    axis.ticks.x = element_blank(),
    axis.ticks.y = element_blank(),
    axis.text.x = element_text(size = 11),
    axis.text.y = element_blank(),
    panel.background = element_blank(),
    panel.grid.major.x = element_line(color = "gray80", size = 0.3),
    legend.position = 'none',
    legend.text = element_text(size=12),
    legend.key= element_blank(),
    plot.margin = margin(20, 20, 20, 20)
  ) +
  coord_cartesian(x = c(0, 3)) +
  facet_wrap(~dist, ncol = 1, scales = "free_y")

Posterior distributions

level_order = c("Pessimistic",
                "Optimistic",
                "Skeptical",
                "Non-informative",
                "Evidence-based"
                )

p1 = samples_posteriors_using_discharge %>%
  
  ggplot(aes(
    # Multiply by 100 to report in percentage
    x = 100 * RD,
    y = factor(`Underlying_Prior`, level = level_order),
    fill = `Underlying_Prior`)) +
  geom_vline(xintercept = seq(-12.5, 12.5, 2.5), linetype = 1, size = 0.3,
             color = "gray80") +
  
  # https://mjskay.github.io/ggdist/articles/slabinterval.html
  stat_halfeye(aes(fill = stat(x < 0)),
               position = "dodge",
               point_interval = median_hdi, # Highest density interval
               .width = c(0.8, 0.95),
               # https://github.com/mjskay/ggdist/issues/70
               interval_size_range = c(1.1,1.9)
  ) +
  scale_fill_manual(values = c("gray85", "#A65041")) +
  labs(x = "\nRisk Difference (%)",
       y = "Underlying Prior\n"
  ) +
  scale_x_continuous(breaks = seq(from = -15, to = 15, 5)) +
  scale_y_discrete(expand = c(0, 0.3)) +
  theme(axis.title.x = element_text(size = 12),
    axis.title.y = element_text(size = 13.5),
    axis.text.y = element_text(size = 12),
    axis.text.x = element_text(size = 12),
    axis.ticks.x = element_blank(),
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    panel.grid.major.x = element_line(color = "gray80", size = 0.3),
    legend.position = "none",
    plot.title.position = 'plot',
    plot.margin = margin(20, 25, 40, 20)
  ) +
  coord_cartesian(xlim = c(-15, 15))

p1_overlay = p1 + inset_element(arrows_plot,
                                ignore_tag = TRUE,
                                align_to = "full",
                                left = unit(4.8, 'cm'),
                                bottom = unit(0, 'cm'),
                                right = unit(16.8, 'cm'),
                                top = unit(1.6, 'cm'))

p1_overlay

Point estimates depict the median. Interval bars depict the 80 and 95% highest density intervals.

# ggsave(width = 7,
#       height = 4.5,
#        here("final_analyses", "output", "plots", "appendix", # File path
#                  "SFigure16_discharge_posteriors_using.pdf")) # File name

t = tibble(
  
  "Underlying Prior" = c(
    "Non-informative",
    "Evidence-based",
    "Skeptical",
    "Optimistic",
    "Pessimistic"
  ),
  
  
  "Pr(< -5%)" = c(round((samples_posteriors_using_discharge %>%
                           filter(`Underlying_Prior` == "Non-informative") %>%
                           summarise(n = sum(RD < -0.05) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_using_discharge %>%
                           filter(`Underlying_Prior` == "Evidence-based") %>%
                           summarise(n = sum(RD < -0.05) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_using_discharge %>%
                           filter(`Underlying_Prior` == "Skeptical") %>%
                           summarise(n = sum(RD < -0.05) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_using_discharge %>%
                           filter(`Underlying_Prior` == "Optimistic") %>%
                           summarise(n = sum(RD < -0.05) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_using_discharge %>%
                           filter(`Underlying_Prior` == "Pessimistic") %>%
                           summarise(n = sum(RD < -0.05) / n()) %>%
                           pull()), 2) * 100
  ),
  
  "Pr(< -2%)" = c(round((samples_posteriors_using_discharge %>%
                           filter(`Underlying_Prior` == "Non-informative") %>%
                           summarise(n = sum(RD < -0.02) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_using_discharge %>%
                           filter(`Underlying_Prior` == "Evidence-based") %>%
                           summarise(n = sum(RD < -0.02) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_using_discharge %>%
                           filter(`Underlying_Prior` == "Skeptical") %>%
                           summarise(n = sum(RD < -0.02) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_using_discharge %>%
                           filter(`Underlying_Prior` == "Optimistic") %>%
                           summarise(n = sum(RD < -0.02) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_using_discharge %>%
                           filter(`Underlying_Prior` == "Pessimistic") %>%
                           summarise(n = sum(RD < -0.02) / n()) %>%
                           pull()), 2) * 100
  ),
  
  "Pr(< -1%)" = c(round((samples_posteriors_using_discharge %>%
                           filter(`Underlying_Prior` == "Non-informative") %>%
                           summarise(n = sum(RD < -0.01) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_using_discharge %>%
                           filter(`Underlying_Prior` == "Evidence-based") %>%
                           summarise(n = sum(RD < -0.01) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_using_discharge %>%
                           filter(`Underlying_Prior` == "Skeptical") %>%
                           summarise(n = sum(RD < -0.01) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_using_discharge %>%
                           filter(`Underlying_Prior` == "Optimistic") %>%
                           summarise(n = sum(RD < -0.01) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_using_discharge %>%
                           filter(`Underlying_Prior` == "Pessimistic") %>%
                           summarise(n = sum(RD < -0.01) / n()) %>%
                           pull()), 2) * 100
  ),
  
  
  "Pr(< 0%)" = c(round((samples_posteriors_using_discharge %>%
                          filter(`Underlying_Prior` == "Non-informative") %>%
                          summarise(n = sum(RD < 0) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_using_discharge %>%
                          filter(`Underlying_Prior` == "Evidence-based") %>%
                          summarise(n = sum(RD < 0) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_using_discharge %>%
                          filter(`Underlying_Prior` == "Skeptical") %>%
                          summarise(n = sum(RD < 0) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_using_discharge %>%
                          filter(`Underlying_Prior` == "Optimistic") %>%
                          summarise(n = sum(RD < 0) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_using_discharge %>%
                          filter(`Underlying_Prior` == "Pessimistic") %>%
                          summarise(n = sum(RD < 0) / n()) %>%
                          pull()), 2) * 100
  ),
  
  "Pr(> 1%)" = c(round((samples_posteriors_using_discharge %>%
                          filter(`Underlying_Prior` == "Non-informative") %>%
                          summarise(n = sum(RD > 0.01) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_using_discharge %>%
                          filter(`Underlying_Prior` == "Evidence-based") %>%
                          summarise(n = sum(RD > 0.01) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_using_discharge %>%
                          filter(`Underlying_Prior` == "Skeptical") %>%
                          summarise(n = sum(RD > 0.01) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_using_discharge %>%
                          filter(`Underlying_Prior` == "Optimistic") %>%
                          summarise(n = sum(RD > 0.01) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_using_discharge %>%
                          filter(`Underlying_Prior` == "Pessimistic") %>%
                          summarise(n = sum(RD > 0.01) / n()) %>%
                          pull()), 2) * 100
  ),
  
  "Pr(> 2%)" = c(round((samples_posteriors_using_discharge %>%
                          filter(`Underlying_Prior` == "Non-informative") %>%
                          summarise(n = sum(RD > 0.02) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_using_discharge %>%
                          filter(`Underlying_Prior` == "Evidence-based") %>%
                          summarise(n = sum(RD > 0.02) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_using_discharge %>%
                          filter(`Underlying_Prior` == "Skeptical") %>%
                          summarise(n = sum(RD > 0.02) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_using_discharge %>%
                          filter(`Underlying_Prior` == "Optimistic") %>%
                          summarise(n = sum(RD > 0.02) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_using_discharge %>%
                          filter(`Underlying_Prior` == "Pessimistic") %>%
                          summarise(n = sum(RD > 0.02) / n()) %>%
                          pull()), 2) * 100
  ),
  
  "Pr(> 5%)" = c(round((samples_posteriors_using_discharge %>%
                          filter(`Underlying_Prior` == "Non-informative") %>%
                          summarise(n = sum(RD > 0.05) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_using_discharge %>%
                          filter(`Underlying_Prior` == "Evidence-based") %>%
                          summarise(n = sum(RD > 0.05) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_using_discharge %>%
                          filter(`Underlying_Prior` == "Skeptical") %>%
                          summarise(n = sum(RD > 0.05) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_using_discharge %>%
                          filter(`Underlying_Prior` == "Optimistic") %>%
                          summarise(n = sum(RD > 0.05) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_using_discharge %>%
                          filter(`Underlying_Prior` == "Pessimistic") %>%
                          summarise(n = sum(RD > 0.05) / n()) %>%
                          pull()), 2) * 100
  )
) 

# https://stackoverflow.com/questions/26548495/reorder-rows-using-custom-order

level_order = c("Evidence-based",
                "Non-informative",
                "Skeptical",
                "Optimistic",
                "Pessimistic"
                )

t = t %>% 
  slice(match(level_order, `Underlying Prior`))

HDI for each posterior distribution

Here are the 95% HDI in risk difference of each distribution.

multiply100 = function(x){x*100}

t1 = samples_posteriors_using_discharge %>%
  select("Underlying_Prior", "RD") %>% 
  # Group by distribution so we can apply median_hdi
  group_by(`Underlying_Prior`) %>% 
  # Calculate the median and 95% quantile interval in log-odds ratio
  ggdist::median_hdi(.width = 0.95) %>% 
  # Select only the relevant columns
  # Arrange order
  arrange(factor(`Underlying_Prior`, levels = c(
    "Evidence-based",
    "Non-informative",
    "Skeptical",
    "Optimistic",
    "Pessimistic"
    ))) %>% 
  select("Underlying_Prior":.upper) %>%
  mutate(across(RD:.upper, ~multiply100(.))) %>%
  mutate(across(RD:.upper, ~round(., 2))) %>%
  rename("Underlying Prior" = Underlying_Prior,
         "Median" = RD,
         "Lower limit" = .lower,
         "Upper limit" = .upper) 

tfinal = left_join(t1, t)

## Joining, by = "Underlying Prior"

tfinal = tfinal %>% 
  gt() %>% 
  tab_spanner(label = "95% HDI",
              columns = c(2:4)) %>% 
  tab_spanner(label = "Propability of Benefit",
              columns = c("Pr(< -5%)", "Pr(< -2%)", "Pr(< -1%)", "Pr(< 0%)")) %>% 
  tab_spanner(label = "Probability of Harm",
              columns = c("Pr(> 1%)", "Pr(> 2%)", "Pr(> 5%)")) %>%
  cols_align(
    align = "center",
    columns = everything()
  ) 
  

tfinal

Underlying Prior	95% HDI			Propability of Benefit				Probability of Harm
	Median	Lower limit	Upper limit	Pr(< -5%)	Pr(< -2%)	Pr(< -1%)	Pr(< 0%)	Pr(> 1%)	Pr(> 2%)	Pr(> 5%)
Evidence-based	-7.87	-11.04	-4.63	96	100	100	100	0	0	0
Non-informative	-8.59	-11.86	-5.30	98	100	100	100	0	0	0
Skeptical	-8.17	-11.33	-4.95	97	100	100	100	0	0	0
Optimistic	-8.43	-11.61	-5.21	98	100	100	100	0	0	0
Pessimistic	-7.89	-11.04	-4.66	96	100	100	100	0	0	0

# tfinal %>% 
#   gtsave(here("final_analyses", "output", "plots", "appendix", # File path
#                   "STable11_discharge_hdi_probabilities_using.png"))

Simple oxygen only

Prior and likelihood distributions

distributions =
  tibble(
    
    "RECOVERY (likelihood)" = list(
      rnorm(
        10e4,
        mean = (posteriors_simple_discharge %>%
                  # Data.mean and data.sd are all the same in every row
                  # So it doesn't matter what row I slice
                  slice(1) %>% 
                  pull(data.mean)),
        sd = (posteriors_simple_discharge %>%
                slice(1) %>% 
                pull(data.sd))
      )),
    
    "Non-informative prior" = list(
      rnorm(
        10e4,
        mean = (posteriors_simple_discharge %>%
                  filter(type == "non-informative") %>%
                  pull(prior.mean)),
        sd = (posteriors_simple_discharge %>%
                filter(type == "non-informative") %>%
                pull(prior.sd))
      )),
    
    "Evidence-based prior" = list(
      rnorm(
        10e4,
        mean = (posteriors_simple_discharge %>%
                  filter(type == "evidence-based") %>%
                  pull(prior.mean)),
        sd = (posteriors_simple_discharge %>%
                filter(type == "evidence-based") %>%
                pull(prior.sd))
      )),
    
    "Skeptical prior" = list(
      rnorm(
        10e4,
        mean = (posteriors_simple_discharge %>%
                  filter(type == "skeptical") %>%
                  pull(prior.mean)),
        sd = (posteriors_simple_discharge %>%
                filter(type == "skeptical") %>%
                pull(prior.sd))
      )),
    
    "Optimistic prior" = list(
      rnorm(
        10e4,
        mean = (posteriors_simple_discharge %>%
                  filter(type == "optimistic") %>%
                  pull(prior.mean)),
        sd = (posteriors_simple_discharge %>%
                filter(type == "optimistic") %>%
                pull(prior.sd))
      )),
    
    "Pessimistic prior" = list(
      rnorm(
        10e4,
        mean = (posteriors_simple_discharge %>%
                  filter(type == "pessimistic") %>%
                  pull(prior.mean)),
        sd = (posteriors_simple_discharge %>%
                filter(type == "pessimistic") %>%
                pull(prior.sd))
      ))
  ) %>% 
  unnest(everything())

# Plot!

distributions %>%
  # make it tidy for ggplot
  pivot_longer(
    cols = c(everything()),
    names_to = "dist", values_to = "samples"
  ) %>%
  # Set order of priors
  mutate(dist = factor(dist, levels =
                         c("RECOVERY (likelihood)",
                           "Non-informative prior",
                           "Evidence-based prior",
                           "Skeptical prior",
                           "Optimistic prior",
                           "Pessimistic prior")
  )
  ) %>%
  # exp() to transform log(OR) to OR
  ggplot(aes(
    x = exp(samples), y = dist,
    fill = dist
  )) +
  
  # https://mjskay.github.io/ggdist/articles/slabinterval.html
  stat_halfeye(position = "dodge",
               point_interval = median_qi, # Quantile interval 
               .width = c(0.95)) +
  geom_vline(xintercept = 1, linetype = 2, size = 0.6) +
  scale_fill_npg()+
  labs(
    x = "\nOdds Ratio",
    y = " "
  ) +
  scale_x_continuous(breaks = seq(from = 0, to = 3, 0.5)) +
  scale_y_discrete(expand = c(0, 0.5)) +
  theme(
    strip.background = element_rect(fill = "#E4E6E7"),
    strip.text.x = element_text(size = 12),
    axis.ticks.x = element_blank(),
    axis.ticks.y = element_blank(),
    axis.text.x = element_text(size = 11),
    axis.text.y = element_blank(),
    panel.background = element_blank(),
    panel.grid.major.x = element_line(color = "gray80", size = 0.3),
    legend.position = 'none',
    legend.text = element_text(size=12),
    legend.key= element_blank(),
    plot.margin = margin(20, 20, 20, 20)
  ) +
  coord_cartesian(x = c(0, 3)) +
  facet_wrap(~dist, ncol = 1, scales = "free_y")

Posterior distributions

level_order = c("Pessimistic",
                "Optimistic",
                "Skeptical",
                "Non-informative",
                "Evidence-based")

p1 = samples_posteriors_simple_discharge %>%
  
  ggplot(aes(
    # Multiply by 100 to report in percentage
    x = 100 * RD,
    y = factor(`Underlying_Prior`, level = level_order),
    fill = `Underlying_Prior`)) +
  geom_vline(xintercept = seq(-12.5, 12.5, 2.5), linetype = 1, size = 0.3,
             color = "gray80") +
  
  # https://mjskay.github.io/ggdist/articles/slabinterval.html
  stat_halfeye(aes(fill = stat(x < 0)),
               position = "dodge",
               point_interval = median_hdi, # Highest density interval
               .width = c(0.8, 0.95),
               # https://github.com/mjskay/ggdist/issues/70
               interval_size_range = c(1.1,1.9)
  ) +
  scale_fill_manual(values = c("gray85", "#9D75B1")) +
  labs(x = "\nRisk Difference (%)",
       y = "Underlying Prior\n"
  ) +
  scale_x_continuous(breaks = seq(from = -15, to = 15, 5)) +
  scale_y_discrete(expand = c(0, 0.3)) +
  theme(axis.title.x = element_text(size = 12),
    axis.title.y = element_text(size = 13.5),
    axis.text.y = element_text(size = 12),
    axis.text.x = element_text(size = 12),
    axis.ticks.x = element_blank(),
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    panel.grid.major.x = element_line(color = "gray80", size = 0.3),
    legend.position = "none",
    plot.title.position = 'plot',
    plot.margin = margin(20, 25, 40, 20)
  ) +
  coord_cartesian(xlim = c(-15, 15))

p1_overlay = p1 + inset_element(arrows_plot,
                                ignore_tag = TRUE,
                                align_to = "full",
                                left = unit(4.8, 'cm'),
                                bottom = unit(0, 'cm'),
                                right = unit(16.8, 'cm'),
                                top = unit(1.6, 'cm'))

p1_overlay

Point estimates depict the median. Interval bars depict the 80 and 95% highest density intervals.

# ggsave(width = 7,
#        height = 4.5,
#         here("final_analyses", "output", "plots", "appendix", # File path
#                   "SFigure17_discharge_posteriors_simple.pdf")) # File name

t = tibble(
  
  "Underlying Prior" = c(
    "Non-informative",
    "Evidence-based",
    "Skeptical",
    "Optimistic",
    "Pessimistic"
  ),
  
  
  "Pr(< -5%)" = c(round((samples_posteriors_simple_discharge %>%
                           filter(`Underlying_Prior` == "Non-informative") %>%
                           summarise(n = sum(RD < -0.05) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_simple_discharge %>%
                           filter(`Underlying_Prior` == "Evidence-based") %>%
                           summarise(n = sum(RD < -0.05) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_simple_discharge %>%
                           filter(`Underlying_Prior` == "Skeptical") %>%
                           summarise(n = sum(RD < -0.05) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_simple_discharge %>%
                           filter(`Underlying_Prior` == "Optimistic") %>%
                           summarise(n = sum(RD < -0.05) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_simple_discharge %>%
                           filter(`Underlying_Prior` == "Pessimistic") %>%
                           summarise(n = sum(RD < -0.05) / n()) %>%
                           pull()), 2) * 100
  ),
  
  "Pr(< -2%)" = c(round((samples_posteriors_simple_discharge %>%
                           filter(`Underlying_Prior` == "Non-informative") %>%
                           summarise(n = sum(RD < -0.02) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_simple_discharge %>%
                           filter(`Underlying_Prior` == "Evidence-based") %>%
                           summarise(n = sum(RD < -0.02) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_simple_discharge %>%
                           filter(`Underlying_Prior` == "Skeptical") %>%
                           summarise(n = sum(RD < -0.02) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_simple_discharge %>%
                           filter(`Underlying_Prior` == "Optimistic") %>%
                           summarise(n = sum(RD < -0.02) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_simple_discharge %>%
                           filter(`Underlying_Prior` == "Pessimistic") %>%
                           summarise(n = sum(RD < -0.02) / n()) %>%
                           pull()), 2) * 100
  ),
  
  "Pr(< -1%)" = c(round((samples_posteriors_simple_discharge %>%
                           filter(`Underlying_Prior` == "Non-informative") %>%
                           summarise(n = sum(RD < -0.01) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_simple_discharge %>%
                           filter(`Underlying_Prior` == "Evidence-based") %>%
                           summarise(n = sum(RD < -0.01) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_simple_discharge %>%
                           filter(`Underlying_Prior` == "Skeptical") %>%
                           summarise(n = sum(RD < -0.01) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_simple_discharge %>%
                           filter(`Underlying_Prior` == "Optimistic") %>%
                           summarise(n = sum(RD < -0.01) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_simple_discharge %>%
                           filter(`Underlying_Prior` == "Pessimistic") %>%
                           summarise(n = sum(RD < -0.01) / n()) %>%
                           pull()), 2) * 100
  ),
  
  
  "Pr(< 0%)" = c(round((samples_posteriors_simple_discharge %>%
                          filter(`Underlying_Prior` == "Non-informative") %>%
                          summarise(n = sum(RD < 0) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_simple_discharge %>%
                          filter(`Underlying_Prior` == "Evidence-based") %>%
                          summarise(n = sum(RD < 0) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_simple_discharge %>%
                          filter(`Underlying_Prior` == "Skeptical") %>%
                          summarise(n = sum(RD < 0) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_simple_discharge %>%
                          filter(`Underlying_Prior` == "Optimistic") %>%
                          summarise(n = sum(RD < 0) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_simple_discharge %>%
                          filter(`Underlying_Prior` == "Pessimistic") %>%
                          summarise(n = sum(RD < 0) / n()) %>%
                          pull()), 2) * 100
  ),
  
  "Pr(> 1%)" = c(round((samples_posteriors_simple_discharge %>%
                          filter(`Underlying_Prior` == "Non-informative") %>%
                          summarise(n = sum(RD > 0.01) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_simple_discharge %>%
                          filter(`Underlying_Prior` == "Evidence-based") %>%
                          summarise(n = sum(RD > 0.01) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_simple_discharge %>%
                          filter(`Underlying_Prior` == "Skeptical") %>%
                          summarise(n = sum(RD > 0.01) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_simple_discharge %>%
                          filter(`Underlying_Prior` == "Optimistic") %>%
                          summarise(n = sum(RD > 0.01) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_simple_discharge %>%
                          filter(`Underlying_Prior` == "Pessimistic") %>%
                          summarise(n = sum(RD > 0.01) / n()) %>%
                          pull()), 2) * 100
  ),
  
  "Pr(> 2%)" = c(round((samples_posteriors_simple_discharge %>%
                          filter(`Underlying_Prior` == "Non-informative") %>%
                          summarise(n = sum(RD > 0.02) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_simple_discharge %>%
                          filter(`Underlying_Prior` == "Evidence-based") %>%
                          summarise(n = sum(RD > 0.02) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_simple_discharge %>%
                          filter(`Underlying_Prior` == "Skeptical") %>%
                          summarise(n = sum(RD > 0.02) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_simple_discharge %>%
                          filter(`Underlying_Prior` == "Optimistic") %>%
                          summarise(n = sum(RD > 0.02) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_simple_discharge %>%
                          filter(`Underlying_Prior` == "Pessimistic") %>%
                          summarise(n = sum(RD > 0.02) / n()) %>%
                          pull()), 2) * 100
  ),
  
  "Pr(> 5%)" = c(round((samples_posteriors_simple_discharge %>%
                          filter(`Underlying_Prior` == "Non-informative") %>%
                          summarise(n = sum(RD > 0.05) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_simple_discharge %>%
                          filter(`Underlying_Prior` == "Evidence-based") %>%
                          summarise(n = sum(RD > 0.05) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_simple_discharge %>%
                          filter(`Underlying_Prior` == "Skeptical") %>%
                          summarise(n = sum(RD > 0.05) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_simple_discharge %>%
                          filter(`Underlying_Prior` == "Optimistic") %>%
                          summarise(n = sum(RD > 0.05) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_simple_discharge %>%
                          filter(`Underlying_Prior` == "Pessimistic") %>%
                          summarise(n = sum(RD > 0.05) / n()) %>%
                          pull()), 2) * 100
  )
) 

# https://stackoverflow.com/questions/26548495/reorder-rows-using-custom-order

level_order = c("Evidence-based",
                "Non-informative",
                "Skeptical",
                "Optimistic",
                "Pessimistic"
                )

t = t %>% 
  slice(match(level_order, `Underlying Prior`))

HDI for each posterior distribution

Here are the 95% HDI in risk difference of each distribution.

multiply100 = function(x){x*100}

t1 = samples_posteriors_simple_discharge %>%
  select("Underlying_Prior", "RD") %>% 
  # Group by distribution so we can apply median_hdi
  group_by(`Underlying_Prior`) %>% 
  # Calculate the median and 95% quantile interval in log-odds ratio
  ggdist::median_hdi(.width = 0.95) %>% 
  # Select only the relevant columns
  # Arrange order
  arrange(factor(`Underlying_Prior`, levels = c(
    "Evidence-based",
    "Non-informative",
    "Skeptical",
    "Optimistic",
    "Pessimistic"
    ))) %>% 
  select("Underlying_Prior":.upper) %>%
  mutate(across(RD:.upper, ~multiply100(.))) %>%
  mutate(across(RD:.upper, ~round(., 2))) %>%
  rename("Underlying Prior" = Underlying_Prior,
         "Median" = RD,
         "Lower limit" = .lower,
         "Upper limit" = .upper) 

tfinal = left_join(t1, t)

## Joining, by = "Underlying Prior"

tfinal = tfinal %>% 
  gt() %>% 
  tab_spanner(label = "95% HDI",
              columns = c(2:4)) %>% 
  tab_spanner(label = "Propability of Benefit",
              columns = c("Pr(< -5%)", "Pr(< -2%)", "Pr(< -1%)", "Pr(< 0%)")) %>% 
  tab_spanner(label = "Probability of Harm",
              columns = c("Pr(> 1%)", "Pr(> 2%)", "Pr(> 5%)")) %>%
  cols_align(
    align = "center",
    columns = everything()
  ) 
  

tfinal

Underlying Prior	95% HDI			Propability of Benefit				Probability of Harm
	Median	Lower limit	Upper limit	Pr(< -5%)	Pr(< -2%)	Pr(< -1%)	Pr(< 0%)	Pr(> 1%)	Pr(> 2%)	Pr(> 5%)
Evidence-based	-6.75	-10.27	-3.08	82	99	100	100	0	0	0
Non-informative	-6.49	-10.25	-2.66	78	99	100	100	0	0	0
Skeptical	-5.92	-9.52	-2.13	68	98	99	100	0	0	0
Optimistic	-6.33	-9.80	-2.53	76	99	100	100	0	0	0
Pessimistic	-5.53	-9.23	-1.78	61	96	99	100	0	0	0

# tfinal %>% 
#   gtsave(here("final_analyses", "output", "plots", "appendix", # File path
#                   "STable12_discharge_hdi_probabilities_simple.png"))

Non-invasive ventilation

Prior and likelihood distributions

distributions =
  tibble(
    
    "RECOVERY (likelihood)" = list(
      rnorm(
        10e4,
        mean = (posteriors_noninvasive_discharge %>%
                  # Data.mean and data.sd are all the same in every row
                  # So it doesn't matter what row I slice
                  slice(1) %>% 
                  pull(data.mean)),
        sd = (posteriors_noninvasive_discharge %>%
                slice(1) %>% 
                pull(data.sd))
      )),
    
    "Non-informative prior" = list(
      rnorm(
        10e4,
        mean = (posteriors_noninvasive_discharge %>%
                  filter(type == "non-informative") %>%
                  pull(prior.mean)),
        sd = (posteriors_noninvasive_discharge %>%
                filter(type == "non-informative") %>%
                pull(prior.sd))
      )),
    
    "Evidence-based prior" = list(
      rnorm(
        10e4,
        mean = (posteriors_noninvasive_discharge %>%
                  filter(type == "evidence-based") %>%
                  pull(prior.mean)),
        sd = (posteriors_noninvasive_discharge %>%
                filter(type == "evidence-based") %>%
                pull(prior.sd))
      )),
    
    "Skeptical prior" = list(
      rnorm(
        10e4,
        mean = (posteriors_noninvasive_discharge %>%
                  filter(type == "skeptical") %>%
                  pull(prior.mean)),
        sd = (posteriors_noninvasive_discharge %>%
                filter(type == "skeptical") %>%
                pull(prior.sd))
      )),
    
    "Optimistic prior" = list(
      rnorm(
        10e4,
        mean = (posteriors_noninvasive_discharge %>%
                  filter(type == "optimistic") %>%
                  pull(prior.mean)),
        sd = (posteriors_noninvasive_discharge %>%
                filter(type == "optimistic") %>%
                pull(prior.sd))
      )),
    
    "Pessimistic prior" = list(
      rnorm(
        10e4,
        mean = (posteriors_noninvasive_discharge %>%
                  filter(type == "pessimistic") %>%
                  pull(prior.mean)),
        sd = (posteriors_noninvasive_discharge %>%
                filter(type == "pessimistic") %>%
                pull(prior.sd))
      ))
  ) %>% 
  unnest(everything())

# Plot!

distributions %>%
  # make it tidy for ggplot
  pivot_longer(
    cols = c(everything()),
    names_to = "dist", values_to = "samples"
  ) %>%
  # Set order of priors
  mutate(dist = factor(dist, levels =
                         c("RECOVERY (likelihood)",
                           "Non-informative prior",
                           "Evidence-based prior",
                           "Skeptical prior",
                           "Optimistic prior",
                           "Pessimistic prior")
  )
  ) %>%
  # exp() to transform log(OR) to OR
  ggplot(aes(
    x = exp(samples), y = dist,
    fill = dist
  )) +
  
  # https://mjskay.github.io/ggdist/articles/slabinterval.html
  stat_halfeye(position = "dodge",
               point_interval = median_qi, # Quantile interval 
               .width = c(0.95)) +
  geom_vline(xintercept = 1, linetype = 2, size = 0.6) +
  scale_fill_npg()+
  labs(
    x = "\nOdds Ratio",
    y = " "
  ) +
  scale_x_continuous(breaks = seq(from = 0, to = 3, 0.5)) +
  scale_y_discrete(expand = c(0, 0.5)) +
  theme(
    strip.background = element_rect(fill = "#E4E6E7"),
    strip.text.x = element_text(size = 12),
    axis.ticks.x = element_blank(),
    axis.ticks.y = element_blank(),
    axis.text.x = element_text(size = 11),
    axis.text.y = element_blank(),
    panel.background = element_blank(),
    panel.grid.major.x = element_line(color = "gray80", size = 0.3),
    legend.position = 'none',
    legend.text = element_text(size=12),
    legend.key= element_blank(),
    plot.margin = margin(20, 20, 20, 20)
  ) +
  coord_cartesian(x = c(0, 3)) +
  facet_wrap(~dist, ncol = 1, scales = "free_y")

Posterior distributions

level_order = c("Pessimistic",
                "Optimistic",
                "Skeptical",
                "Non-informative",
                "Evidence-based")

p1 = samples_posteriors_noninvasive_discharge %>%
  
  ggplot(aes(
    # Multiply by 100 to report in percentage
    x = 100 * RD,
    y = factor(`Underlying_Prior`, level = level_order),
    fill = `Underlying_Prior`)) +
  geom_vline(xintercept = seq(-12.5, 12.5, 2.5), linetype = 1, size = 0.3,
             color = "gray80") +
  
  # https://mjskay.github.io/ggdist/articles/slabinterval.html
  stat_halfeye(aes(fill = stat(x < 0)),
               position = "dodge",
               point_interval = median_hdi, # Highest density interval
               .width = c(0.8, 0.95),
               # https://github.com/mjskay/ggdist/issues/70
               interval_size_range = c(1.1,1.9)
  ) +
  scale_fill_manual(values = c("gray85", "#CE9F5B")) +
  labs(x = "\nRisk Difference (%)",
       y = "Underlying Prior\n"
  ) +
  scale_x_continuous(breaks = seq(from = -15, to = 15, 5)) +
  scale_y_discrete(expand = c(0, 0.3)) +
  theme(axis.title.x = element_text(size = 12),
    axis.title.y = element_text(size = 13.5),
    axis.text.y = element_text(size = 12),
    axis.text.x = element_text(size = 12),
    axis.ticks.x = element_blank(),
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    panel.grid.major.x = element_line(color = "gray80", size = 0.3),
    legend.position = "none",
    plot.title.position = 'plot',
    plot.margin = margin(20, 25, 40, 20)
  ) +
  coord_cartesian(xlim = c(-15, 15))

p1_overlay = p1 + inset_element(arrows_plot,
                                ignore_tag = TRUE,
                                align_to = "full",
                                left = unit(4.8, 'cm'),
                                bottom = unit(0, 'cm'),
                                right = unit(16.8, 'cm'),
                                top = unit(1.6, 'cm'))

p1_overlay

Point estimates depict the median. Interval bars depict the 80 and 95% highest density intervals.

# ggsave(width = 7,
#        height = 4.5,
#         here("final_analyses", "output", "plots", "appendix", # File path
#                   "SFigure18_discharge_posteriors_noninvasive.pdf")) # File name

t = tibble(
  
  "Underlying Prior" = c(
    "Non-informative",
    "Evidence-based",
    "Skeptical",
    "Optimistic",
    "Pessimistic"
  ),
  
  
  "Pr(< -5%)" = c(round((samples_posteriors_noninvasive_discharge %>%
                           filter(`Underlying_Prior` == "Non-informative") %>%
                           summarise(n = sum(RD < -0.05) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_noninvasive_discharge %>%
                           filter(`Underlying_Prior` == "Evidence-based") %>%
                           summarise(n = sum(RD < -0.05) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_noninvasive_discharge %>%
                           filter(`Underlying_Prior` == "Skeptical") %>%
                           summarise(n = sum(RD < -0.05) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_noninvasive_discharge %>%
                           filter(`Underlying_Prior` == "Optimistic") %>%
                           summarise(n = sum(RD < -0.05) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_noninvasive_discharge %>%
                           filter(`Underlying_Prior` == "Pessimistic") %>%
                           summarise(n = sum(RD < -0.05) / n()) %>%
                           pull()), 2) * 100
  ),
  
  "Pr(< -2%)" = c(round((samples_posteriors_noninvasive_discharge %>%
                           filter(`Underlying_Prior` == "Non-informative") %>%
                           summarise(n = sum(RD < -0.02) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_noninvasive_discharge %>%
                           filter(`Underlying_Prior` == "Evidence-based") %>%
                           summarise(n = sum(RD < -0.02) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_noninvasive_discharge %>%
                           filter(`Underlying_Prior` == "Skeptical") %>%
                           summarise(n = sum(RD < -0.02) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_noninvasive_discharge %>%
                           filter(`Underlying_Prior` == "Optimistic") %>%
                           summarise(n = sum(RD < -0.02) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_noninvasive_discharge %>%
                           filter(`Underlying_Prior` == "Pessimistic") %>%
                           summarise(n = sum(RD < -0.02) / n()) %>%
                           pull()), 2) * 100
  ),
  
  "Pr(< -1%)" = c(round((samples_posteriors_noninvasive_discharge %>%
                           filter(`Underlying_Prior` == "Non-informative") %>%
                           summarise(n = sum(RD < -0.01) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_noninvasive_discharge %>%
                           filter(`Underlying_Prior` == "Evidence-based") %>%
                           summarise(n = sum(RD < -0.01) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_noninvasive_discharge %>%
                           filter(`Underlying_Prior` == "Skeptical") %>%
                           summarise(n = sum(RD < -0.01) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_noninvasive_discharge %>%
                           filter(`Underlying_Prior` == "Optimistic") %>%
                           summarise(n = sum(RD < -0.01) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_noninvasive_discharge %>%
                           filter(`Underlying_Prior` == "Pessimistic") %>%
                           summarise(n = sum(RD < -0.01) / n()) %>%
                           pull()), 2) * 100
  ),
  
  
  "Pr(< 0%)" = c(round((samples_posteriors_noninvasive_discharge %>%
                          filter(`Underlying_Prior` == "Non-informative") %>%
                          summarise(n = sum(RD < 0) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_noninvasive_discharge %>%
                          filter(`Underlying_Prior` == "Evidence-based") %>%
                          summarise(n = sum(RD < 0) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_noninvasive_discharge %>%
                          filter(`Underlying_Prior` == "Skeptical") %>%
                          summarise(n = sum(RD < 0) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_noninvasive_discharge %>%
                          filter(`Underlying_Prior` == "Optimistic") %>%
                          summarise(n = sum(RD < 0) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_noninvasive_discharge %>%
                          filter(`Underlying_Prior` == "Pessimistic") %>%
                          summarise(n = sum(RD < 0) / n()) %>%
                          pull()), 2) * 100
  ),
  
  "Pr(> 1%)" = c(round((samples_posteriors_noninvasive_discharge %>%
                          filter(`Underlying_Prior` == "Non-informative") %>%
                          summarise(n = sum(RD > 0.01) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_noninvasive_discharge %>%
                          filter(`Underlying_Prior` == "Evidence-based") %>%
                          summarise(n = sum(RD > 0.01) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_noninvasive_discharge %>%
                          filter(`Underlying_Prior` == "Skeptical") %>%
                          summarise(n = sum(RD > 0.01) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_noninvasive_discharge %>%
                          filter(`Underlying_Prior` == "Optimistic") %>%
                          summarise(n = sum(RD > 0.01) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_noninvasive_discharge %>%
                          filter(`Underlying_Prior` == "Pessimistic") %>%
                          summarise(n = sum(RD > 0.01) / n()) %>%
                          pull()), 2) * 100
  ),
  
  "Pr(> 2%)" = c(round((samples_posteriors_noninvasive_discharge %>%
                          filter(`Underlying_Prior` == "Non-informative") %>%
                          summarise(n = sum(RD > 0.02) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_noninvasive_discharge %>%
                          filter(`Underlying_Prior` == "Evidence-based") %>%
                          summarise(n = sum(RD > 0.02) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_noninvasive_discharge %>%
                          filter(`Underlying_Prior` == "Skeptical") %>%
                          summarise(n = sum(RD > 0.02) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_noninvasive_discharge %>%
                          filter(`Underlying_Prior` == "Optimistic") %>%
                          summarise(n = sum(RD > 0.02) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_noninvasive_discharge %>%
                          filter(`Underlying_Prior` == "Pessimistic") %>%
                          summarise(n = sum(RD > 0.02) / n()) %>%
                          pull()), 2) * 100
  ),
  
  "Pr(> 5%)" = c(round((samples_posteriors_noninvasive_discharge %>%
                          filter(`Underlying_Prior` == "Non-informative") %>%
                          summarise(n = sum(RD > 0.05) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_noninvasive_discharge %>%
                          filter(`Underlying_Prior` == "Evidence-based") %>%
                          summarise(n = sum(RD > 0.05) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_noninvasive_discharge %>%
                          filter(`Underlying_Prior` == "Skeptical") %>%
                          summarise(n = sum(RD > 0.05) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_noninvasive_discharge %>%
                          filter(`Underlying_Prior` == "Optimistic") %>%
                          summarise(n = sum(RD > 0.05) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_noninvasive_discharge %>%
                          filter(`Underlying_Prior` == "Pessimistic") %>%
                          summarise(n = sum(RD > 0.05) / n()) %>%
                          pull()), 2) * 100
  )
) 

# https://stackoverflow.com/questions/26548495/reorder-rows-using-custom-order

level_order = c("Evidence-based",
                "Non-informative",
                "Skeptical",
                "Optimistic",
                "Pessimistic"
                )

t = t %>% 
  slice(match(level_order, `Underlying Prior`))

HDI for each posterior distribution

Here are the 95% HDI in risk difference of each distribution.

multiply100 = function(x){x*100}

t1 = samples_posteriors_noninvasive_discharge %>%
  select("Underlying_Prior", "RD") %>% 
  # Group by distribution so we can apply median_hdi
  group_by(`Underlying_Prior`) %>% 
  # Calculate the median and 95% quantile interval in log-odds ratio
  ggdist::median_hdi(.width = 0.95) %>% 
  # Select only the relevant columns
  # Arrange order
  arrange(factor(`Underlying_Prior`, levels = c(
    "Evidence-based",
    "Non-informative",
    "Skeptical",
    "Optimistic",
    "Pessimistic"
    ))) %>% 
  select("Underlying_Prior":.upper) %>%
  mutate(across(RD:.upper, ~multiply100(.))) %>%
  mutate(across(RD:.upper, ~round(., 2))) %>%
  rename("Underlying Prior" = Underlying_Prior,
         "Median" = RD,
         "Lower limit" = .lower,
         "Upper limit" = .upper) 

tfinal = left_join(t1, t)

## Joining, by = "Underlying Prior"

tfinal = tfinal %>% 
  gt() %>% 
  tab_spanner(label = "95% HDI",
              columns = c(2:4)) %>% 
  tab_spanner(label = "Propability of Benefit",
              columns = c("Pr(< -5%)", "Pr(< -2%)", "Pr(< -1%)", "Pr(< 0%)")) %>% 
  tab_spanner(label = "Probability of Harm",
              columns = c("Pr(> 1%)", "Pr(> 2%)", "Pr(> 5%)")) %>%
  cols_align(
    align = "center",
    columns = everything()
  ) 
  

tfinal

Underlying Prior	95% HDI			Propability of Benefit				Probability of Harm
	Median	Lower limit	Upper limit	Pr(< -5%)	Pr(< -2%)	Pr(< -1%)	Pr(< 0%)	Pr(> 1%)	Pr(> 2%)	Pr(> 5%)
Evidence-based	-7.13	-11.74	-2.56	82	99	100	100	0	0	0
Non-informative	-7.21	-12.09	-2.49	82	98	99	100	0	0	0
Skeptical	-6.61	-11.22	-2.05	75	98	99	100	0	0	0
Optimistic	-7.06	-11.65	-2.47	81	98	100	100	0	0	0
Pessimistic	-6.14	-10.71	-1.53	69	96	99	100	0	0	0

# tfinal %>% 
#   gtsave(here("final_analyses", "output", "plots", "appendix", # File path
#                   "STable13_discharge_hdi_probabilities_noninvasive.png"))

Invasive mechanical ventilation

Prior and likelihood distributions

distributions =
  tibble(
    
    "RECOVERY (likelihood)" = list(
      rnorm(
        10e4,
        mean = (posteriors_invasive_discharge %>%
                  # Data.mean and data.sd are all the same in every row
                  # So it doesn't matter what row I slice
                  slice(1) %>% 
                  pull(data.mean)),
        sd = (posteriors_invasive_discharge %>%
                slice(1) %>% 
                pull(data.sd))
      )),
    
    "Non-informative prior" = list(
      rnorm(
        10e4,
        mean = (posteriors_invasive_discharge %>%
                  filter(type == "non-informative") %>%
                  pull(prior.mean)),
        sd = (posteriors_invasive_discharge %>%
                filter(type == "non-informative") %>%
                pull(prior.sd))
      )),
    
    "Evidence-based prior" = list(
      rnorm(
        10e4,
        mean = (posteriors_invasive_discharge %>%
                  filter(type == "evidence-based") %>%
                  pull(prior.mean)),
        sd = (posteriors_invasive_discharge %>%
                filter(type == "evidence-based") %>%
                pull(prior.sd))
      )),
    
    "Skeptical prior" = list(
      rnorm(
        10e4,
        mean = (posteriors_invasive_discharge %>%
                  filter(type == "skeptical") %>%
                  pull(prior.mean)),
        sd = (posteriors_invasive_discharge %>%
                filter(type == "skeptical") %>%
                pull(prior.sd))
      )),
    
    "Optimistic prior" = list(
      rnorm(
        10e4,
        mean = (posteriors_invasive_discharge %>%
                  filter(type == "optimistic") %>%
                  pull(prior.mean)),
        sd = (posteriors_invasive_discharge %>%
                filter(type == "optimistic") %>%
                pull(prior.sd))
      )),
    
    "Pessimistic prior" = list(
      rnorm(
        10e4,
        mean = (posteriors_invasive_discharge %>%
                  filter(type == "pessimistic") %>%
                  pull(prior.mean)),
        sd = (posteriors_invasive_discharge %>%
                filter(type == "pessimistic") %>%
                pull(prior.sd))
      ))
  ) %>% 
  unnest(everything())

# Plot!

distributions %>%
  # make it tidy for ggplot
  pivot_longer(
    cols = c(everything()),
    names_to = "dist", values_to = "samples"
  ) %>%
  # Set order of priors
  mutate(dist = factor(dist, levels =
                         c("RECOVERY (likelihood)",
                           "Non-informative prior",
                           "Evidence-based prior",
                           "Skeptical prior",
                           "Optimistic prior",
                           "Pessimistic prior")
  )
  ) %>%
  # exp() to transform log(OR) to OR
  ggplot(aes(
    x = exp(samples), y = dist,
    fill = dist
  )) +
  
  # https://mjskay.github.io/ggdist/articles/slabinterval.html
  stat_halfeye(position = "dodge",
               point_interval = median_qi, # Quantile interval 
               .width = c(0.95)) +
  geom_vline(xintercept = 1, linetype = 2, size = 0.6) +
 scale_fill_npg()+
  labs(
    x = "\nOdds Ratio",
    y = " "
  ) +
  scale_x_continuous(breaks = seq(from = 0, to = 3, 0.5)) +
  scale_y_discrete(expand = c(0, 0.5)) +
  theme(
    strip.background = element_rect(fill = "#E4E6E7"),
    strip.text.x = element_text(size = 12),
    axis.ticks.x = element_blank(),
    axis.ticks.y = element_blank(),
    axis.text.x = element_text(size = 11),
    axis.text.y = element_blank(),
    panel.background = element_blank(),
    panel.grid.major.x = element_line(color = "gray80", size = 0.3),
    legend.position = 'none',
    legend.text = element_text(size=12),
    legend.key= element_blank(),
    plot.margin = margin(20, 20, 20, 20)
  ) +
  coord_cartesian(x = c(0, 3)) +
  facet_wrap(~dist, ncol = 1, scales = "free_y")

Posterior distributions

level_order = c("Pessimistic",
                "Optimistic",
                "Skeptical",
                "Non-informative",
                "Evidence-based")

p1 = samples_posteriors_invasive_discharge %>%
  
  ggplot(aes(
    # Multiply by 100 to report in percentage
    x = 100 * RD,
    y = factor(`Underlying_Prior`, level = level_order),
    fill = `Underlying_Prior`)) +
  geom_vline(xintercept = seq(-12.5, 12.5, 2.5), linetype = 1, size = 0.3,
             color = "gray80") +
  
  # https://mjskay.github.io/ggdist/articles/slabinterval.html
  stat_halfeye(aes(fill = stat(x < 0)),
               position = "dodge",
               point_interval = median_hdi, # Highest density interval
               .width = c(0.8, 0.95),
               # https://github.com/mjskay/ggdist/issues/70
               interval_size_range = c(1.1,1.9)
  ) +
  scale_fill_manual(values = c("gray85", "#5891BF")) +
  labs(x = "\nRisk Difference (%)",
       y = "Underlying Prior\n"
  ) +
  scale_x_continuous(breaks = seq(from = -15, to = 15, 5)) +
  scale_y_discrete(expand = c(0, 0.3)) +
  theme(axis.title.x = element_text(size = 12),
    axis.title.y = element_text(size = 13.5),
    axis.text.y = element_text(size = 12),
    axis.text.x = element_text(size = 12),
    axis.ticks.x = element_blank(),
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    panel.grid.major.x = element_line(color = "gray80", size = 0.3),
    legend.position = "none",
    plot.title.position = 'plot',
    plot.margin = margin(20, 25, 40, 20)
  ) +
  coord_cartesian(xlim = c(-15, 15))

p1_overlay = p1 + inset_element(arrows_plot,
                                ignore_tag = TRUE,
                                align_to = "full",
                                left = unit(4.8, 'cm'),
                                bottom = unit(0, 'cm'),
                                right = unit(16.8, 'cm'),
                                top = unit(1.6, 'cm'))

p1_overlay

Point estimates depict the median. Interval bars depict the 80 and 95% highest density intervals.

# ggsave(width = 7,
#        height = 4.5,
#         here("final_analyses", "output", "plots", "appendix", # File path
#                   "SFigure19_discharge_posteriors_invasive.pdf")) # File name

t = tibble(
  
  "Underlying Prior" = c(
    "Non-informative",
    "Evidence-based",
    "Skeptical",
    "Optimistic",
    "Pessimistic"
  ),
  
  
  "Pr(< -5%)" = c(round((samples_posteriors_invasive_discharge %>%
                           filter(`Underlying_Prior` == "Non-informative") %>%
                           summarise(n = sum(RD < -0.05) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_invasive_discharge %>%
                           filter(`Underlying_Prior` == "Evidence-based") %>%
                           summarise(n = sum(RD < -0.05) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_invasive_discharge %>%
                           filter(`Underlying_Prior` == "Skeptical") %>%
                           summarise(n = sum(RD < -0.05) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_invasive_discharge %>%
                           filter(`Underlying_Prior` == "Optimistic") %>%
                           summarise(n = sum(RD < -0.05) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_invasive_discharge %>%
                           filter(`Underlying_Prior` == "Pessimistic") %>%
                           summarise(n = sum(RD < -0.05) / n()) %>%
                           pull()), 2) * 100
  ),
  
  "Pr(< -2%)" = c(round((samples_posteriors_invasive_discharge %>%
                           filter(`Underlying_Prior` == "Non-informative") %>%
                           summarise(n = sum(RD < -0.02) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_invasive_discharge %>%
                           filter(`Underlying_Prior` == "Evidence-based") %>%
                           summarise(n = sum(RD < -0.02) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_invasive_discharge %>%
                           filter(`Underlying_Prior` == "Skeptical") %>%
                           summarise(n = sum(RD < -0.02) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_invasive_discharge %>%
                           filter(`Underlying_Prior` == "Optimistic") %>%
                           summarise(n = sum(RD < -0.02) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_invasive_discharge %>%
                           filter(`Underlying_Prior` == "Pessimistic") %>%
                           summarise(n = sum(RD < -0.02) / n()) %>%
                           pull()), 2) * 100
  ),
  
  "Pr(< -1%)" = c(round((samples_posteriors_invasive_discharge %>%
                           filter(`Underlying_Prior` == "Non-informative") %>%
                           summarise(n = sum(RD < -0.01) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_invasive_discharge %>%
                           filter(`Underlying_Prior` == "Evidence-based") %>%
                           summarise(n = sum(RD < -0.01) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_invasive_discharge %>%
                           filter(`Underlying_Prior` == "Skeptical") %>%
                           summarise(n = sum(RD < -0.01) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_invasive_discharge %>%
                           filter(`Underlying_Prior` == "Optimistic") %>%
                           summarise(n = sum(RD < -0.01) / n()) %>%
                           pull()), 2) * 100,
                  
                  round((samples_posteriors_invasive_discharge %>%
                           filter(`Underlying_Prior` == "Pessimistic") %>%
                           summarise(n = sum(RD < -0.01) / n()) %>%
                           pull()), 2) * 100
  ),
  
  
  "Pr(< 0%)" = c(round((samples_posteriors_invasive_discharge %>%
                          filter(`Underlying_Prior` == "Non-informative") %>%
                          summarise(n = sum(RD < 0) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_invasive_discharge %>%
                          filter(`Underlying_Prior` == "Evidence-based") %>%
                          summarise(n = sum(RD < 0) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_invasive_discharge %>%
                          filter(`Underlying_Prior` == "Skeptical") %>%
                          summarise(n = sum(RD < 0) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_invasive_discharge %>%
                          filter(`Underlying_Prior` == "Optimistic") %>%
                          summarise(n = sum(RD < 0) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_invasive_discharge %>%
                          filter(`Underlying_Prior` == "Pessimistic") %>%
                          summarise(n = sum(RD < 0) / n()) %>%
                          pull()), 2) * 100
  ),
  
  "Pr(> 1%)" = c(round((samples_posteriors_invasive_discharge %>%
                          filter(`Underlying_Prior` == "Non-informative") %>%
                          summarise(n = sum(RD > 0.01) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_invasive_discharge %>%
                          filter(`Underlying_Prior` == "Evidence-based") %>%
                          summarise(n = sum(RD > 0.01) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_invasive_discharge %>%
                          filter(`Underlying_Prior` == "Skeptical") %>%
                          summarise(n = sum(RD > 0.01) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_invasive_discharge %>%
                          filter(`Underlying_Prior` == "Optimistic") %>%
                          summarise(n = sum(RD > 0.01) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_invasive_discharge %>%
                          filter(`Underlying_Prior` == "Pessimistic") %>%
                          summarise(n = sum(RD > 0.01) / n()) %>%
                          pull()), 2) * 100
  ),
  
  "Pr(> 2%)" = c(round((samples_posteriors_invasive_discharge %>%
                          filter(`Underlying_Prior` == "Non-informative") %>%
                          summarise(n = sum(RD > 0.02) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_invasive_discharge %>%
                          filter(`Underlying_Prior` == "Evidence-based") %>%
                          summarise(n = sum(RD > 0.02) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_invasive_discharge %>%
                          filter(`Underlying_Prior` == "Skeptical") %>%
                          summarise(n = sum(RD > 0.02) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_invasive_discharge %>%
                          filter(`Underlying_Prior` == "Optimistic") %>%
                          summarise(n = sum(RD > 0.02) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_invasive_discharge %>%
                          filter(`Underlying_Prior` == "Pessimistic") %>%
                          summarise(n = sum(RD > 0.02) / n()) %>%
                          pull()), 2) * 100
  ),
  
  "Pr(> 5%)" = c(round((samples_posteriors_invasive_discharge %>%
                          filter(`Underlying_Prior` == "Non-informative") %>%
                          summarise(n = sum(RD > 0.05) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_invasive_discharge %>%
                          filter(`Underlying_Prior` == "Evidence-based") %>%
                          summarise(n = sum(RD > 0.05) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_invasive_discharge %>%
                          filter(`Underlying_Prior` == "Skeptical") %>%
                          summarise(n = sum(RD > 0.05) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_invasive_discharge %>%
                          filter(`Underlying_Prior` == "Optimistic") %>%
                          summarise(n = sum(RD > 0.05) / n()) %>%
                          pull()), 2) * 100,
                 
                 round((samples_posteriors_invasive_discharge %>%
                          filter(`Underlying_Prior` == "Pessimistic") %>%
                          summarise(n = sum(RD > 0.05) / n()) %>%
                          pull()), 2) * 100
  )
) 

# https://stackoverflow.com/questions/26548495/reorder-rows-using-custom-order

level_order = c("Evidence-based",
                "Non-informative",
                "Skeptical",
                "Optimistic",
                "Pessimistic"
                )

t = t %>% 
  slice(match(level_order, `Underlying Prior`))

HDI for each posterior distribution

Here are the 95% HDI in risk difference of each distribution.

multiply100 = function(x){x*100}

t1 = samples_posteriors_invasive_discharge %>%
  select("Underlying_Prior", "RD") %>% 
  # Group by distribution so we can apply median_hdi
  group_by(`Underlying_Prior`) %>% 
  # Calculate the median and 95% quantile interval in log-odds ratio
  ggdist::median_hdi(.width = 0.95) %>% 
  # Select only the relevant columns
  # Arrange order
  arrange(factor(`Underlying_Prior`, levels = c(
    "Evidence-based",
    "Non-informative",
    "Skeptical",
    "Optimistic",
    "Pessimistic"
    ))) %>% 
  select("Underlying_Prior":.upper) %>%
  mutate(across(RD:.upper, ~multiply100(.))) %>%
  mutate(across(RD:.upper, ~round(., 2))) %>%
  rename("Underlying Prior" = Underlying_Prior,
         "Median" = RD,
         "Lower limit" = .lower,
         "Upper limit" = .upper) 

tfinal = left_join(t1, t)

## Joining, by = "Underlying Prior"

tfinal = tfinal %>% 
  gt() %>% 
  tab_spanner(label = "95% HDI",
              columns = c(2:4)) %>% 
  tab_spanner(label = "Propability of Benefit",
              columns = c("Pr(< -5%)", "Pr(< -2%)", "Pr(< -1%)", "Pr(< 0%)")) %>% 
  tab_spanner(label = "Probability of Harm",
              columns = c("Pr(> 1%)", "Pr(> 2%)", "Pr(> 5%)")) %>%
  cols_align(
    align = "center",
    columns = everything()
  ) 
  

tfinal

Underlying Prior	95% HDI			Propability of Benefit				Probability of Harm
	Median	Lower limit	Upper limit	Pr(< -5%)	Pr(< -2%)	Pr(< -1%)	Pr(< 0%)	Pr(> 1%)	Pr(> 2%)	Pr(> 5%)
Evidence-based	-3.96	-10.33	1.58	37	75	85	92	4	1	0
Non-informative	-3.43	-10.64	2.87	33	66	77	85	8	4	0
Skeptical	-2.48	-8.27	2.93	20	57	70	82	9	4	0
Optimistic	-3.36	-9.53	2.15	30	68	80	89	5	2	0
Pessimistic	-1.64	-7.28	3.54	13	45	59	73	15	8	0

# tfinal %>% 
#   gtsave(here("final_analyses", "output", "plots", "appendix", # File path
#                   "STable14_discharge_hdi_probabilities_invasive.png"))

sessionInfo()

## R version 4.0.4 (2021-02-15)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] Cairo_1.5-12.2   ggsci_2.9        patchwork_1.1.1  ggdist_2.4.0    
##  [5] janitor_2.1.0    metafor_2.4-0    Matrix_1.3-2     here_1.0.1      
##  [9] flextable_0.6.5  gt_0.3.0         kableExtra_1.3.4 forcats_0.5.1   
## [13] stringr_1.4.0    dplyr_1.0.5      purrr_0.3.4      readr_1.4.0     
## [17] tidyr_1.1.3      tibble_3.1.1     ggplot2_3.3.3    tidyverse_1.3.1 
## [21] rio_0.5.26       magrittr_2.0.1   plyr_1.8.6       pacman_0.5.1    
## 
## loaded via a namespace (and not attached):
##  [1] nlme_3.1-152          fs_1.5.0              lubridate_1.7.10     
##  [4] webshot_0.5.2         httr_1.4.2            rprojroot_2.0.2      
##  [7] tools_4.0.4           backports_1.2.1       bslib_0.2.4          
## [10] utf8_1.2.1            R6_2.5.0              DBI_1.1.1            
## [13] colorspace_2.0-1      withr_2.4.2           tidyselect_1.1.1     
## [16] curl_4.3.1            compiler_4.0.4        cli_2.5.0            
## [19] rvest_1.0.0           HDInterval_0.2.2      xml2_1.3.2           
## [22] officer_0.3.18        labeling_0.4.2        sass_0.3.1           
## [25] checkmate_2.0.0       scales_1.1.1          systemfonts_1.0.1    
## [28] digest_0.6.27         foreign_0.8-81        rmarkdown_2.8        
## [31] svglite_2.0.0         base64enc_0.1-3       pkgconfig_2.0.3      
## [34] htmltools_0.5.1.1     highr_0.9             dbplyr_2.1.1         
## [37] fastmap_1.1.0         rlang_0.4.11          readxl_1.3.1         
## [40] rstudioapi_0.13       farver_2.1.0          jquerylib_0.1.4      
## [43] generics_0.1.0        jsonlite_1.7.2        distributional_0.2.2 
## [46] zip_2.1.1             Rcpp_1.0.6            munsell_0.5.0        
## [49] fansi_0.4.2           gdtools_0.2.3         lifecycle_1.0.0      
## [52] stringi_1.6.1         yaml_2.2.1            snakecase_0.11.0     
## [55] grid_4.0.4            crayon_1.4.1          lattice_0.20-41      
## [58] haven_2.4.1           hms_1.0.0             knitr_1.33           
## [61] pillar_1.6.0          uuid_0.1-4            reprex_2.0.0         
## [64] glue_1.4.2            evaluate_0.14         V8_3.4.0             
## [67] data.table_1.14.0     modelr_0.1.8          vctrs_0.3.8          
## [70] cellranger_1.1.0      gtable_0.3.0          assertthat_0.2.1     
## [73] cachem_1.0.4          xfun_0.22             openxlsx_4.2.3       
## [76] broom_0.7.6           viridisLite_0.4.0     memoise_2.0.0        
## [79] ellipsis_0.3.2        rdocsyntax_0.4.1.9000